A Density Functional Study of Tautomerism of Uracil and Cytosine

Apr 15, 1994 - Received: November 17, 1993; In Final Form: March 2, 1994'. The structure .... vog Molecule ... TABLE 2 Relative Energies (kJ/mol) of t...
2 downloads 0 Views 951KB Size
Download PDF
J. Phys. Chem. 1994,98, 5653-5660

5653

A Density Functional Study of Tautomerism of Uracil and Cytosine Dario A. Estrin, Luca Paglieri, and Giorgina Corongiu' Centro di Ricerca, Sviluppo e Studi Superiori in Sardegna (CRS4),P.O.Box 488, 09100 Cagliari, Italy Received: November 17, 1993; In Final Form: March 2, 1994'

The structure and the relative energies of all possible tautomeric forms of the uracil and cytosine molecules have been determined using both local and gradient-corrected density functional methods. The calculations have been performed with double-zeta plus polarization basis sets and the geometries optimized with analytic gradient techniques. The vibrational frequencies and the contribution of the zero-point energies have also been computed. In the uracil case, the dioxo form is predicted to be the most stable. In the cytosine case, three tautomers are found to be very close in energy with the oxo-amino form slightly more stable. The infrared absorption intensities and frequencies for the uracil and the two more stable tautomeric forms of the cytosine molecules are reported and compared with experimental spectra. The agreement with experiment and correlated ab-initio methods is good for geometries, energetics, and vibrational frequencies.

1. Introduction

stabilities of the different tautomers of the uracil and cytosine molecules. Since in both cases some of the tautomers are very A large amount of work has been performed on the tautomerism close in energy, this study can represent an additional test for the of nucleic acid bases using both theoreticalld and e ~ p e r i m e n t a l ~ - ~ DFT methods applicability. In order to assess the relative stability approaches. Much of the interest is due to the fact that tautomers of the tautomers, we have computed also the zero-point energy induce alterations in the normal base pairing, leading to the contribution to the internal energy. Vibrational frequencies and possibility of spontaneous mutationsl0J1 in the DNA or R N A infrared absorption intensities are compared with experimental helices. IR spectra. The theoretical determination of tautomeric stabilities is achieved by evaluating differences between electronic energies, 2. Computational Methodology including the difference in the zero-point vibrational energy. This is numerically a very demanding task, since it involves differences The calculations in this study have been carried out with the Molecole-DFT program,Igwhere the Kohn-Sham self-consistent between two large and almost equal quantities. From a procedureZois applied for obtaining the electronic density and methodology viewpoint, it has been shown that the evaluation of relative stabilities for tautomers requires an adequate treatment energy through the determination of a set of one-electron orbitals of correlation effects*,6needed to properly describe the alternate (the solution of the Kohn-Sham equations). Gaussian basis sets multiple bond structure of these molecules. The use of correlated are used for the expansion of the one-electron orbitals. The methods has been usually limited to single-point calculations on electronic density is also expanded in an additional Gaussian basis sets2' The coefficients for the fit of the electronic density structures optimized with lower levels of t h e ~ r y . l - ~It is well are obtained by minimizing the error in the Coulomb repulsion known, however, that the Hartree-Fock determined geometries energy. The use of this fit reduces the integral evaluation process show considerable deviations from (gas phase) experimental dependency from N4to WM (whereNis the number of functions values, and therefore it is necessary to perform the structure optimization with post-Hartree-Fock methods in order to obtain in the orbital basis set and M is the total number of functions in a more satisfactory agreement.I* Since traditional ab-initio the auxiliary basis, typically of comparable size to N). methods which include correlation effects scale as Nm with m > Matrix elements of the exchange-correlation potential are 5 ( N is the number of basis functions of the system), the evaluated by a numerical integration scheme based on the grids corresponding calculations are computationally very demanding and quadratures reported by Becke.*2 In performing the numeriand impose serious limitations on the size of the system. cal integration of the exchange-correlation potential matrix elements, in the S C F step, we have used a grid with 25 radial An alternative approach to conventional ab-initio methods is shells for the carbon, nitrogen, and oxygen atoms, 20 for the based on the use of density functional theory (DFT), which has hydrogen atoms, and 50 angular points per shell. At the end of been applied extensively to the study of solids and surfaces13J4 the S C F procedure, the exchange-correlation energy and its and only recently is used in chemical problem~.l5-~~ The density contribution to the gradients have been evaluated using an functional methods are very interesting from the computational augmented grid, with 35 radial shells for the carbon, nitrogen, point of view, since the algorithm scales as N or less in some and oxygen atoms, 30 for the hydrogen atoms, and 194 angular types of implementations, and it accounts for electron correlation points per radial shell. effects. It has also been shownI6 that DFT methods can handle problems that, in the ab-initio formalism, should be treated with The double-zeta plus polarization orbital basis sets optimized multiconfigurational approaches, since it takes into account most by Sim et al.23~2~ for LDA-DFT calculations have been used. The of the dynamical and nondynamical correlation effects, still contraction pattern is (521 1/41 1/1) for carbon nitrogen, and keeping the simplicity of the independent particle model. This oxygen and (41/1) for hydrogen. This corresponds, as quality, is important in the case of the uracil and cytosine tautomerism, to the 6-31G** basis set used in standard ab-initio methods. For since the presence of alternatedouble bonds can make inadequate the electronic density expansion basis set the contraction pattern those approximations which include only dynamical correlation. was (1 111111/ 1 1/ 1) and (1 11111/ 1) for the carbon, nitrogen, oxygen, and hydrogen atoms. In this work we will demonstrate the reliability of the density functional approach in predicting the the structure and relative An extended description of the technical detail of the program is given in ref 19. The restricted S C F procedure has been applied in all cases, and no symmetry has been imposed on the systems. 0 Abstract published in Advance ACS Abstracts, April 15, 1994. 0022-3654f 94f 2098-5653$04.50f 0

0 1994 American Chemical Society

5654

Estrin et al.

The Journal of Physical Chemistry, Vol. 98, No. 22, 1994

9

Y

O8

O8

vog rr

u2

UI

U3

H -08

C6

u4

us

U6

Figure 1. Six tautomers of uracil.

The computations have been performed at two different levels of the density functional theory. In one case we have used the local density approximation (LDA), with the Gunnarsson and Lundqvist25 parametrization of the homogeneous electron gas exchange-correlation energy density. In the second case two different gradient-corrected density functionals have been considered. In both functionals the correlation term is given by the VoskoZ6 parametrization of the homogeneous electron gas with the Perdew2' gradient corrections. Regarding the gradient corrections in the exchange term, in one case the PerdewZ8 expression has been used, while in the other case the one proposed by BeckeZ9 has been used. The results obtained will be labeled LDA, NL1, and NL2. The geometry optimization has been performed by using a quasi-Newton minimization method in Cartesian coordinates with analytical energy gradients. The Cartesian force constants and dipole moment derivatives are calculated by numerical differentiation of the energy gradients and dipole moments using Cartesian displacements of 0.015 au. From the force constant matrix, vibrational frequencies and zero-point energies have been obtained. The infrared absorption intensities have been calculated by taking the numerical derivatives of the dipole moment and transforming them to the corresponding ones with respect to the normal m0des.3~In order to evaluate the relative energies of the different tautomers, the internal energy of a particular molecule has been evaluated by adding to the total electronic energy the zero-point energy.

3. Results 3.1. Uracil Tautomers. The uracil molecule may exist in various tautomeric forms differing from each other by the position of the protons, which may be bound to either ring nitrogen atoms or oxygen atoms. The six uracil tautomers object of the present study are shown in Figure 1. The optimized geometries have been obtained for all tautomers using both the local (LDA) and the two different gradient-corrected density functionals ( N L 1 and NL2) mentioned in the previous section. We report in Table 1 only the structural parameters corresponding to the more stable dioxo form, U1,since it is the only one for which experimental information is available. In Table 1 our results are compared with those obtained with the HF-offset forces optimization method: a HF-3-21G3' geometry optimization, and the X-ray crystallographic results.32 Even if intermolecular hydrogen bonds may produce a noticeable distortion in the molecular ge0metry,~3 nevertheless bond distances and angles obtained from X-rays experiments may still be used as a guide. The results from Table 1 show that the different density functionals provide very similar geometries which compare favorably with the other theoretical and the experimental results.

TABLE 1: Structural Parameters (A and deg) for the Uracil Molecule HFC exptd parameter LDA' NL14 NL2' HFb r(N 1C2) r(c2N3j r(N3C4) r(C4C5) ric5~6j r(C6N1) r(C207) r(C408) r(N1H) r(N3H) r(C5H) r(C6H)

1.382 1.371 1.398 1.443 1.347 1.358 1.216 1.221 1.026 1.030 1.096 1.098

1.408 1.396 1.428 1.465 1.363 1.383 1.230 1.235 1.028 1.031 1.098 1.100

1.404 1.392 1.422 1.463 1.361 1.378 1.228 1.233 1.027 1.030 1.097 1.100

1.383 1.378 1.403 1.457 1.348 1.382 1.214 1.215 1.004 1.007 1.076 1.079

1.380 1.374 1.396 1.457 1.326 1.375 1.211 1.212 0.998 1.002 1.066 1.069

1.374 1.381 1.380 1.444 1.343 1.370 1.219 1.233

L(NlC2N3) L(C2N3C4) L(N3C4C5) L(C4C5C6) L(C5C6Nl) L(C6NlC2) L(NlC207) L(N3C408) L(C6NlH) L(C2N3H) L(C4C5H) L(C5C6H)

113.2 127.8 114.2 119.2 122.2 123.5 122.5 120.0 121.4 116.1 118.0 122.2

112.9 127.7 114.1 119.4 122.3 123.5 122.6 120.0 121.0 116.1 118.3 122.4

113.0 127.8 114.0 119.4 122.2 123.6 122.6 120.1 121.0 116.0 118.2 122.3

113.7 127.6 114.3 119.2 122.0 123.1 122.6 120.1 121.3 115.8 118.5 122.8

113.3 128.1 113.5 119.8 122.4 122.9 122.9 120.7 120.8 115.6 118.4 121.9

115.4 126.4 114.1 120.7 121.2 122.0 122.9 120.5

This work. Reference 6. Results obtained with a 6-31G* basis set. Reference 3 1. Resultsobtained with a 3-21G basis set. Reference 32.

TABLE 2 Relative Energies (kJ/mol) of the Uracil Tautomers; Comparison with Other Theoretical and Experimental Results method U1 U2 U3 U4 u5 LDA' LDA ZPE' NLl' NLl + ZPE' ~ ~ 2 4 NL2 ZPE' MP2b MP3b MP4(SDQ)b MP4(SDQ)C CISDC expt

+

+

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

45.3 44.1 46.9 45.6 46.3 45.1 46.0 44.1 54.7 50.4 56.8

46.3 45.2 49.2 47.9 48.5 47.2 52.3 51.6 61.5

50.1 49.0 55.4 54.0 53.7 52.3 50.4 47.7 66.9 59.6 71.5 79 f 2Sd

82.1 80.0 83.0 80.6 82.9 80.4 82.9 77.7 88.4

U6 87.4 84.2 86.1 82.9 85.6 82.3 97.2 93.2 104.3

92 f 42d

This work. Reference6. Results obtainedwith a 6-3 1G**basis set at HF/6-31G* geometries. Reference 2. Results obtained with a 6-31G** basis set at HF/6-31G** geometries. Reference 8. (I

However, the DFT methods overestimate the bond distances involving hydrogen atoms by about 0.02 A; this has already been pointed out for other LDA calc~lations,3~ and it is probably due to inaccuracies in the density functionals used; the use of gradientcorrected density functionals does not seem to improve the situation. By comparing the DFT results with those of refs 6 and 32 for the other five tautomers, we find the same kind of trends: good agreement for the bond distances between atoms in the ring and an overestimation of bond distances involving hydrogen atoms. As is expected, density functional optimized geometries show less bond localization than traditional ab-initio methods such as H F and MP2. The geometrical information related to the U2U6 tautomers is made available as supplementary material (Tables IS-VS). In Table 2 we report the tautomers' energies relative to the lowest energy value of tautomer U 1 (dioxo form) with and without the zero-point energy corrections (ZPE). In agreement with previous calculations$ in all cases the zero-point energy contribution destabilizes the dioxo form, and its magnitude, even if not at all negligible, does not change the predicted ordering. Our estimated zero-point corrections range from 1.0 to 3.5 kJ/mol, in disagreement with previous results,2 obtained with a HF 6-3 1G** calculation, which yield zero-point energy differences

Tautomerism of Uracil and Cytosine

The Journal of Physical Chemistry, Vol. 98, No. 22, 1994 5655

TABLE 3: Dipole Moment for the Tautomers of Uracil; Comparison with Other Theoretical and Experimental R d t s method U1 U2 U3 U4 US U6 LDA" 4.51 3.35 4.98 1.25 6.71 7.39 NL" 4.43 3.20 4.92 1.29 6.54 7.12 NL24 4.44 3.23 4.93 1.27 6.56 1.16 HFb H Fc MPZd

4.8 4.72 4.21 4.73 3.86 4.16

3.4 3.49

5.2 5.25

1.6 1.31

6.5

7.6

LDA' exptf exptg This work. Reference 3 1. Results obtained with a 3-21G basis set. Reference 2. HF/6-31G** at HF/6-31G** geometries. Reference 2. MP2/6-31GZ at MP2/6-31G* geometry. Reference 36. Local density approximation;Hedin-Lundqvist exchange-correlation potential. Double numerical + polarization basis set. 1 Reference 37. g Reference 38.

not larger than 0.3 kJ/mol. The predicted range in the zeropoint energy differences is not very much affected by the level of DFT used. We found the same ordering ( U l , U2, U3, U4, US, U6) for both the local and the two gradient-corrected density functionals, with very similar relative energies. Even the LDA, known to yield unreliable results, for example, in the energy changes associated with chemical reactions or bond dissociations,34 does a very good job in the evaluation of relative tautomeric energies. All previous theoretical calculations reported the dioxo tautomer U1 as the most stable. However, the relative stability of the other tautomers is very sensitive to both the basis set and the correlation energy estimation. Tautomers U3 and U4 are very close in energy, and for example the tautomer U4 is predicted to be the more stable in HF, MP2, and MP3 calculations6 with a 6-3 1G** basis set but less stable in an MP4 calculation.6 These are theresults60fsingle-point computations witha 6-31G** basis set carried out on geometries optimized at the H F level with a 6-31G* basis set. At the H F level, using a 3-21G basis set, the U3 tautomer is more stable than U4, contrary to the HF6-3 1G** results.31 Our results are in agreement with the most accurate calculation so far reported6 (MP4-6-31G**) regarding to the tautomers energy ordering, but in general, our computed relative energies are smaller. For example, MP4 results6 predict that the tautomer U3 is 6.8 kJ/mol less stable than tautomer U2, while in the DFT studies this difference is 1.l, 2.2, and 1.9 kJ/mol for LDA, N L l , and NL2 density functionals, respectively. Also, the DFT calculations place the tautomer U6 much closer in energy to tautomer US than the MP4 results. The experimental values* are available only for the tautomers U3 and U4. Since the associated errors are very large, these results are far from conclusive, and further investigation on the subject is necessary in order to draw definitive conclusions. The prediction of accurate dipole moments is a very important issue, because the magnitude of the dipole moment is strongly related to the tautomeric stability in polar environments. In Table 3 we gathered the values of the calculated dipole moments for the six uracil tautomers and the available experimental values. The calculated dipole moments of all tautomers show very little sensitivity to the density functional used, with the results of the two gradient-corrected density functionals almost identical. The useof additional setsof polarization functions in DFTcalculations has shown to improve the agreement with the experiment35 for the water and the water dimer systems; however, in this case a calculation using three sets of d polarization functions for the C, 0,and N atoms and three sets of p polarization functions for the hydrogen atoms produced very little changes in the calculated dipole moment (less than 0.1 D). It is expected, according to the results obtained, that the ordering of tautomers U2 and U3 will

be reversed in water solution, because of the higher dipole moment of tautomer U3 and the very small magnitude of the energy difference U2-U3 in gas phase. Vibrational frequencies and integrated absorption intensities have been calculated for all uracil tautomers, but we report only the results corresponding to the most stable tautomer, since it is unlikely that the higher-energy tautomers can be detected experimentally. The information regarding tautomers U2, U3, U4, U5, and U6 is made available as supplementary material (Tables VIS-XS). The three density functionals results and the experimental value^^^.^^ for the U1 tautomer are reported in Table 4. We do not include the HF results of ref 41, since they were obtained by scaling with an empirically adjusted parameter. The normalmodes assignment has been performed by graphically displaying each normal mode using the program GMC.42 The results obtained by using the LDA and the gradient-corrected density functionals have comparable errors with respect to the experimental results, with average percentage deviations of the calculated frequencies of 2.7,3.2, and 3.2 for the LDA, NL1, and NL2 density functionals, respectively, when using the experimental data of ref 40. If we make the comparison with the experimental frequencies of ref 39, then the average percentage deviations are 2.9, 3.0, and 2.9, respectively. Our assignments of the vibrational normal modes are in agreement with those of Chin et al.39b and Could et a1.,41 with the exception of modes u9, vIo, and uI I . For these three modes we observe some discrepancies among the LDA and the gradient-corrected density functionals. In addition, discrepancies exist among other published re~ults;3~~40 therefore, it is difficult to make a definitive assignment of these peaks. It must be pointed out that anharmonic contributions are present in the experimental frequencies, and probably the excellent agreement found between computed harmonic frequencies and experimental frequencies is to some degree due to a cancellation of errors. The accurate prediction of infrared absorption intensities is, as well as its experimental measurement, a very difficult task, and the reported resultsmust betaken withcautionsincevariations are observed when taking the spectra in different matrices. In order to judge the general accuracy of the calculation, a simulated spectrum (using the NL2 density functional) is reported in Figure 2, which was calculated using a Lorentian function with a constant half-width of 10 cm-1 to represent each peak, with the peak height given by the computed intensity. Also, a simulated spectrum constructed from the experimentally determined frequencies and relative intensities of ref 39 is shown for comparison. 3.2. Cytosine Tautomers. Cytosine may exist in different tautomeric forms, according to the position of a proton, which may bind either to a ring nitrogen atom or to the exocyclic oxygen or amine nitrogen atoms (Figure 3). The first three tautomers (Cl, C2, and C3 in Figure 3) have been extensively st~diedl,~-3~.43+ because they are very close energetically, and the stability order has proved to be extremely sensitive to the level of theory used. Experimentally, it has been found that cytosine exists as a mixture of the amino-oxo and amino-hydroxy tautomers, C1 and C3, respectively, in an Ar or Nz low-temperature matrix.43 The imine tautomer, C2, has not been identified, even if it is also very close energetically to the parent amino-oxo compound. The structural parameters of the cytosine molecule, C1, are reported in Table 5 and compared with the results of a H F 3-21G c a l c ~ l a t i o n a, ~ ~ HF and MP2 6-31G* calc~lation,4~ and X-ray experimental data.32 With the exception of ref 45, all calculations performed on cytosine tautomerism assumed all tautomers to be planar. We found the hydrogen atoms of the amino group in tautomers C1, C3,andC4 to beout oftheplane. Allplanarstructurescorrespond to transition states, very close in energy to the nonplanar minimum-

5656 The Journal of Physical Chemistry, Vol. 98, No. 22, 1994

Estrin et al.

TABLE 4 Calculated Frequencies (cm-l) and IR Intensities (km mol-') for the Uracil Molecule; Comparison with Experimenhl Frequencies and Relative Intensities mode In-Plane ModesC V I N1H s v 2 N3H s v 3 C5H s v 4 C6H s V J C 2 0 S, N3H N l H b

LDA

NLl

NL2

exDt'

eXDtb ~

C 4 0 S, N3H C5H C6H b C5C6 C 4 0 S, N 1H b N l H b, N1C6 C 2 0 s v g N1H C 2 0 N3H C 4 0 b v10 C5H C6H b in ph N3H b, rings v12 ring s, C5H b vi3 C5H C6H b, ooph, N l H bring s V I 4 C5H b, ring s v i 5 ring def V I 6 ring def V I 7 ring breathing v18 ring def v i 9 ring def, C 2 0 C 4 0 b v20 ring def, C 4 0 b V ~ CI 4 0 C 2 0 b v7

vg

Out-of-Plane Modesc ~ 2 C6H 2 C5H b ooph v23 C 4 0 b, C5H C6H b in ph ~ 2C 4 20 b V ~ C5H J C6H b, ring def V26 N3H b, ring def V27 N l H b, ring def v28 ring def v29 ring def V30 ring def

3582 (138) 3539 (74) 3227 (6) 3184 (0) 1813 (705)

3546 (98) 3507 (46) 3196 (1) 3153 (3) 1720 (724)

3563 (105) 3526 (50) 3206 (1) 3157 (2) 1745 (713)

3485 3435

1768 (669)

1677 (637)

1702 (636)

1706

1678 (30) 1521 (107) 1405 (17) 1424 (81) 1357 (1.4) 1260 (14) 1179 (69) 1083 (8) 994 (6) 960 (6) 778 (2) 545 (5) 529 (2) 505 (19) 316 (22)

1618 (8) 1477 (41) 1403 (6) 1345 (93) 1367 (66) 1210 (6) 1154 (102) 1058 (5) 973 (7) 929 (14) 745 (3) 533 (5) 519 (5) 501 (21) 367 (18)

1631 (19) 1486 (50) 1408 (5) 1375 (94) 1352 (60) 1215 (2) 1164 (104) 1066 (6) 974 (7) 934 (12) 749 (3) 535 (5) 519 (5) 500 (20) 368 (19)

1643 1472 1400 1389 1359 1217 1185 1075 980 958 759 562 537 516 39 1

926 (1) 786 (66) 747 (37) 697 (22) 663 (67) 580 (32) 388 (20) 160 (0) 137 (2)

928 (1) 775 (50) 707 (33) 687 (15) 648 (78) 552 (30) 379 (21) 155 (0) 128 (1)

921 (1) 772 (52) 710 (35) 689 (1 5) 646 (79) 552 (29) 377 (22) 153 (0) 130 (1)

987 804 757 718 662 55 1 41 1 185

~~~

3482 (166) 3433 (100) 3130 (4) 2970 (8) 1774d(680) 1762d 1758d 1733d 1720"(291) 1707d 1699d 1644 (33) 1473 (83) 1461 (56) 1401 (56) 1389 (21) 1219 (4) 1186 (109) 1076 (14) 963 (2) 958 (7) 719 (12) 557 (17) 537 (7) 516 (23) 393 (33)

1764

806 (175) 769 (125) 685 (14) 664 (100) 557 (25)

*

Reference 40. Reference 39. b = bending, s = stretching, in ph = in phase, ooph = out of phase, def = deformation. d Frequencies refer todifferent observed maxima; intensity refers to the whole Uracil (Ul) Infrared Spectrum 15

DFT.NL2 vs. Exp.

I

"'W "I0

I

- DFT-NL2 ExP 07

07 I O

L

i

I

CI

I'WH'O 8: 00

i

500

loo0

1500

2MM

25M)

3000

..

3500

c4 4000

Frequency (a")

Figure 2. Uracil (Ul) infrared spectrum;comparisonof DFT-NL2 results with experimental values.

energy structures. In the case of the cytosine molecule, we found a barrier of 0.3 and 1.3 kJ/mol using local and nonlocal density functionals, respectively. Computations carried out on the vibrational frequencies yield only real frequencies for the nonplanar structures, whereas one imaginary frequency, corresponding to the pyramidalization of the amino groups, is found for the planar structures. Our DFT calculations show a good agreement with both experiment and theoretical results. As pointed out previously,

"a: c3

c2

""ij N8

H

07

NI

07

H @ :

07

c5

07

C6

Figure 3. Six tautomers of cytosine.

all three density functionals provide approximately the same geometries. Hydrogen bonds are, as usual, overestimated by about 0.02 A. As in the uracil case, the DFToptimized structures show less localization of bonds than the HF optimized structures. This can be particularly important, for example, in the study of the amino-hydroxy tautomer, C3, for which two resonant structures, of approximately equal weight, exist for this molecule, and using a H F based method, only one of them is obtained. The optimized structures of the C2, C3, C4, C5, and C6 tautomers are made available as supplementary material (Tables XIS-XIVS).

Tautomerism

of Uracil and Cytosine

The Journal of Physical Chemistry, Vol. 98, No. 22, 1994 5657

TABLE 5: Structural Parameters (A and deg) for the Cytosine Molecule

TABLE 7: Dipole Moment for the Tautomers of Cytosine; Comparison with Other Theoretical and Experimental Results

LDA'

NLl'

NL2'

HFb

HFc MPZC exptld

method

c1

C2

C3

C4

C5

C6

r(N 1C2) r(C2N3) r(N3C4) r(C4C5) r(C5C6) r(C6N1) r(C207) r(C4N8) r(N 1H) r(N8H9) r(N8HlO) r(C5H) r(C6H)

1.414 1.358 1.315 1.426 1.356 1.342 1.221 1.353 1.026 1.023 1.018 1.097 1.098

1.443 1.382 1.339 1.447 1.371 1.364 1.235 1.383 1.027 1.025 1.022 1.100 1.101

1.438 1.380 1.330 1.444 1.369 1.360 1.232 1.379 1.027 1.025 1.021 1.098 1.099

1.415 1.369 1.298 1.443 1.337 1.354 1.211 1.344 0.998 0.998 0.995 1.067 1.070

1.402 1.362 1.296 1.444 1.339 1.348 1.196 1.347 0.996 0.995 0.993 1.071 1.074

1.418 1.377 1.319 1.437 1.357 1.356 1.226 1.357 1.014 1.010 1.007 1.083 1.086

1.392 1.358 1.339 1.433 1.357 1.360 1.237 1.324

LDA' NL' NL2' HFb HFc MP2C LDAd exptf

6.71 6.51 6.51 7.2 7.08 7.40 6.66,6.93 7.0

4.83 4.74 4.74 5.3

3.51 3.43 3.42 3.7

8.49 8.13 8.14 8.6

1.78 1.72 1.73 1.8

5.17 5.01 5.06 4.5

L(NlC2N3) L(C2N3C4) L(N3C4C5) L(C4C5C6) L(C5C6N 1) L(C6NlC2) L(NlC207) L(N3C4N8) L(C6N 1H) L(C4N8H9) ~ ( c 4 N 8 H10) L(C4C5H) L(C5C6H)

116.2 120.6 123.7 115.9 120.1 123.5 118.0 117.6 121.7 117.5 120.3 121.0 123.5

116.3 120.1 124.2 116.2 120.0 123.2 118.4 116.7 121.5 118.0 121.1 121.5 122.9

116.3 120.2 124.1 116.1 120.2 123.1 118.2 116.9 121.6 118.0 121.2 121.9 122.6

115.2 122.1 122.8 116.1 120.8 123.0 118.9 118.3 121.3 117.9 122.5 121.7 122.4

116.4 120.4 124.0 115.6 120.3 123.3 118.3 117.6

116.0 119.9 124.5 116.0 119.6 124.0 125.3 116.7

118.6 120.5 121.5 117.0 121.2 121.2 118.9 118.3

117.4 121.1 122.2 123.1

117.9 122.5 122.5 123.4

parameter

~~

+

L(C5C4N8H91 174.5 169.5 169.3 180.0 172.6 L(CSC4N8HlO) 10.2 20.9 21.2 0.0 13.9 a This work. Reference 31. Planar geometry; results obtained with a 3-21G basis set. CReference 45. All results calculated with a 6-31G* basis set. The MP2 results have been obtained assuming a planar geometry. Reference 32.

TABLE 6 Relative Energies (kJ/mol) of Six Cytosine Tautomers; Comparison with Other Theoretical and Experimental Results method

C1

C2

C3

C4

C5

C6

LDA' LDA + ZPE' NLl' NL1 ZPE' NL2' NL2 + ZPE' HFb H Fc MP2c CCd expt'

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

7.5 10.3 5.7 7.6 6.4 8.1 1.7 2.7 5.8 2.2

1.5 2.2 4.0 4.1 2.9 3.2 15.9 2.5 1.0 -4.2

30.1 28.7 26.6 25.2 26.5 25.4 29.8

55.5 57.4 55.3 56.2 55.4 56.2 76.2

78.8 79.8 76.8 76.4 77.3 77.0 94.8

+

a This work. Reference 31. Results obtained with a 3-21G basis set. CReference 45. Results obtained with a 6-31G* basis set. Planar geometry. Reference 36. Local density approximation; Hedin-Lundqvist exchange-correlationpotential. The first value represents the double numerical (d] and the second, the double numerical polarization basis set. f Reference 38.

=O

a This work. Reference 3 1. Planar geometries;results obtained with a 3-21G basis set. Reference 1. Results obtained with a 6-31G* basis set at 3-21G geometries. Reference 44. Coupled cluster method with SDT excitations using first-order correlation orbitals. DZP basis set. Results include ZPE calculated at HF/3-21G level and scaled by a 0.91 factor. Reference 43.

The relative stability of the six cytosine tautomers is presented

in Table 6, with and without the inclusion of the zero-point energy contributions. In contrast to the uracil tautomerism, the inclusion of zero-point energies can be very important since the C1, C2, and C3 tautomers a r e very close energetically, and small changes may substantially alter the predicted order. The order we obtain, without including zero-point energies, is C1, C3, C2, C4, C5, C 6 for all three density functionals considered. The order of tautomers C 3 and C 2 is opposite to the one predicted by the HF 3-21G~alculation3~ but agrees with6-31G* HFandMP2results.I The same order was found also from computations with the coupled cluster method (with a DZP basis set).44 However, in ref 44 the tautomer C 3 is reported to be slightly more stable than C1, even a t the HF level, whereas a HF 6-31G* calculation predicted the C1 tautomer to be more stable. According to ref 44, C 3 is found

+

to be more stable than C 1 if larger basis sets are used with inclusion of p polarization functions on the hydrogen atoms. The DFT computations predict the C 1 tautomer to be themost stabledespite the use of a reasonable basis set with polarization functions. For the higher-energy tautomers, no calculations with correlation effects have so far been reported. Regarding the contributions of the zero-point energy, several calculations have been reported lately. Could and Hillier46 obtained with a HF 6-31G** calculation the value of 0.9 kJ/mol for the C3-C1 tautomer difference, while Les and Adamowicz,s with a HF 3-21G calculation, obtained a value of -2.5 kJ/mol. These authors claim that a value of -1.5 kJ/mol can be inferred from the IR experimental information given by Szczesniak et aL43 Since several frequencies were not observed in the experimental spectra, because of the very low intensity of the bands, the only way of gathering zero-point energies from experimental spectra is in using calculated frequencies when experimental information is not available. This procedure, however, can give only a qualitative rather than quantitative estimate, since it depends upon the level of theory used for the calculation. For the C3-C1 tautomer difference we obtain a value of 0.6,0.1, and 0.3 kJ/mol using the LDA, the NLI, and the N L 2 density functionals, respectively. The results are closer to those obtained with HF 6-31G** calculations46 than those with HF 3-21Gcalculations. Since thedifferencesarevery small, small variations in the geometry may cause the discrepancy, and the assumption of a planar geometry in the 3-21G calculation may affect the results. The difference C2-CI in our case is found to be positive, destabilizing the C2 form, in agreement with the previous calculations. In all the other cases, zero-point energy differences a r e very small, and their contributions do not affect the ordering found. Experimentally, a t 15 K i n argon and nitrogen mat rice^,"^ both tautomers C 1 and C3 coexist, implying an almost zero free energy difference. The C3 is found to be slightly more abundant, whereas the C2 tautomer is not found a t all. Our results agree with the experiments, but we predict the tautomer C1 to be the most stable. Taking as 0.0 the experimental C3-C1 energy difference, we have in the three density functionals considered an error smaller than 1 kcal/mol, which is comparable to the highest level ab-initio calculation reported.44 In Table 7 we gather the values of the calculated and experimental dipole moments for the cytosine tautomers. The results obtained did not change substantially by performing the calculations with three sets of polarization functions, instead of one, as for the uracil case. Our results agree reasonably well with experiment in the case of C1. It is expected than in water solution the tautomer C 4 could be observed, since its high dipole moment guarantees a large solvent stabilization. Vibrational frequencies as well as infrared absorption intensities have been calculated for all cytosine tautomers, since the vibrational analysis has been used for evaluating the nuclear motion contributions to the tautomers stability. Results regarding

5658

The Journal of Physical Chemistry, Vol. 98, No. 22, 1994

Estrin et al.

TABLE 8: Calculated Frequencies (cm-l) and JR Interlsities (km mol-’) for Amino-Oxo Tautomer of Cytosine Molecule; Comparison with Experimental Frequencies and Relative Intensities mode’ LDA NLl NL2 exptb 3682 (70) 3608 (38) 3621 (41) 3565 (81) V I N8H9 N8H10 s ooph NlH s N8H9 N8H10 s in ph u4 C5H s u5 C6H s U6 c 2 0 s u7 C5C6 N3C4 s vg N8H2 sciss, C4N8 s v9 N3C4 C4C5 S, N l H b uI0ring s, N4H9 N4H10 C5H C6H b U I IN3C4 C4C5 S, N1H b v12 C5H C6H bin ph, N l H b ~ 1 N3C2 3 C4N8 s ~ 1 C5H 4 C6H b ooph, N l H b ~ 1 C5C6 5 C6N1 s V I 6 ring wag in pl, CO b Ul7 ring def in pl v l g NlC2 C4C5 s, ring def ~ 1 C5H 9 C6H b ooph 00pl v2

v3

ring breathing b 00p1 b 00pl u23 C5H N l H b oopl, ring def oopl u20 ~ 2C 1 20 N l H ~ 2 C4N8 2 C5H ~ 2N 4 1H

v25 V26 v27

v2g V29

b 00pl

ring s ring s N8H2 tors, C20 b C 2 0 C4N8 b

ring def oopl

V ~ N8H2 O

C207 b

3551 (114) 3527 (64) 3190 (0) 3171 (1) 1779 (816) 1704 (388) 1601 (308) 1589 (140) 1509 (26) 1440 (44) 1349 (30) 1302 (34) 1186 (36) 1107 (4) 1081 (36) 989 (1) 921 (0) 926 (1) 778 (3) 766 (31) 732 (2) 696 (48) 633 (49) 563 (2) 532 (3) 525 (8) 513 (6) 396 (13) 337 (4) 166 (140) 193 (62) 123 (3)

3519 (67) 3477 (37) 3149 (3) 3134 (4) 1686 (755) 1642 (390) 1601 (117) 1527 (99) 1470 (154) 1412 (63) 1325 (75) 1221 (18) 1184 (45) 1097 (4) 1070 (59) 970 (1) 919 (0) 898 (3) 741 ( 5 ) 735 (29) 722 (8) 689 (35) 615 ( 5 5 ) 557 (2) 528 ( 5 ) 512 (6) 500 (40) 385 (23) 348 (100) 337 (112) 192 (8) 123 (3)

3529 (73) 3487 (39) 3165 (3) 3148 (3) 1709 (752) 1653 (367) 1599 (146) 1538 (111) 1479 (139) 1417 (64) 1330 (65) 1228 (20) 1187 (45) 1101 (4) 1074 (57) 971 (1) 918 (0) 902 (3) 744 ( 5 ) 737 (28) 723 (11) 692 (33) 611 (59) 556 (2) 527 (6) 513 (9) 502 (35) 380 (23) 356 (184) 339 (27) 191 (8) 120 (2)

3471 (143) 3441 (209) 1720 (872) 1656 (424) 1595 (71) 1539 (97) 1475 (189) 1422 (56) 1337 (46) 1244 (36) 1192 (70) 1124 (36) 1090 (53) 1087 (8) 818 (14) 781 (41) 747 ( 5 ) 637 (8) 614 (40) 623 (6) 575 (6) 537 (6) 532 ( 5 ) 409 (20) 360 (8) 330 (30) 232e 197c

b oopl b def oopl a b = bending, s = stretching, wag = wagging, sciss = scissoring, def = deformation, tors = torsion, in ph = in phase, ooph = out of phase, oopl = out of ring plane. b Reference 43. Frequencies for crystalline solid. u31 C 2 0 C4N8 ~ 3 C4N8 2 00p1 v33 ring

the infrared spectra are reported only for the two more stable tautomers, C 1 and C3. Information regarding other tautomers is made available as supplementary material (Tables XVSXVIIIS). The results corresponding to the cytosine molecule, C1, are reported in Table 8. All three density functionals considered give an excellent agreement with experimental re~ults.~3The average percentage deviations in the calculated frequencies from the experimental results are 4.9, 3.6, and 3.3 for the LDA, the NL1, and the NL2 density functionals, respectively. The results corresponding to the C3 tautomer are given in Table 9. The average percentage deviations in this case are 3.9, 2.5, and 2.5 for the LDA, NL1, and NL2 density functionals, respectively. As expected, even if the LDA yields satisfactory results, the gradient-corrected density functionals provide better agreement with experiments. The assignments of the vibrational modes are in general agreement with those reported in ref 43; some discrepancies are found in some low-frequencey out-of-plane modes, which are known to be very difficult to assign. Additional discrepancies may be due to the fact that the D I T geometries are only approximately planar for the ring atoms and the H of the amino group is out of the plane. On the other hand, planar structures have been found to be saddle points, having an imaginary frequency corresponding to the pyramidalization of the amino group for both the cytosine and the amino-hydroxy tautomers, C1 and C3, whereas the assignments in ref 43 were based on a H F 3-21G scaled calculation, assuming planar geometries for C1 and C3. To judge the general accuracy of the calculation, simulated spectra for C 1 and C3 tautomers are reported in Figures 4 and 5, respectively, and compared with those constructed from the relative intensities and frequencies reported in ref 43. The spectra have been constructed as in the uracil case. Since all three density

functional results yield similar agreement with the experimental spectra, only the results obtained with the NL2 functional are reported. An excellent overall agreement is obtained, without the necessity of any scaling procedure.

4. Conclusions From the examples discussed in this paper, it appears that current implementations of density functional theory are capable to describe the phenomena of tautomerism reasonably well. Both local and nonlocal functionals yield results which are comparable to those obtained by using correlated ab-initio techniques. The LDA gives satisfactory structural properties and vibrational frequencies. The lack of conclusive experimental or high level ab-initio information regarding the relative stabilities of the tautomers makes it difficult to determine which type of density functional gives the most accurate answers. Ab-initiocalculations including dynamical as well as nondynamical (near-degeneracy) correlation effects are needed to “calibrate” the reliability of the density functional theory in treating this kind of problem. However, it seems that, unlike other chemical problems, for example bond dissociation energies, the use of nonlocal corrections does not improve substantially the quality of the results obtained. Indeed, it can be considered as encouraging the fact that the three density functionals used in this work yield very close answers: both gradient-corrected density functionals gave almost identical results regarding geometries and very similar results concerning relative energies and vibrational frequencies. The DFT methodology is highly competitive because of its computational efficiency. To provide data on the computational effort associated with the calculations presented in this work, in Table 10 we report CPU times, on an IBM RS6000-560 workstation, for the uracil and cytosine molecules (without symmetry and using the grids and basis sets described in section

Tautomerism of Uracil and Cytosine

The Journal of Physical Chemistry, Vol. 98, No. 22, 1994 5659

TABLE 9 Calculated Frequencies (cm-’) and IR Intensities (km mol-’) for Amino-Hydroxy Tautomer of Cytosine Molecule; Comparison with Experimental Frequencies and Relative Intensities modea

LDA

NL1

NL2

exptb

3687 (69) 3647 (100) 3538 (73) 3169 (3) 3171 (1) 1682 (614) 1608 (124) 1642 (258) 1529 (24) 1491 (308) 1422 (101) 1328 (33) 1385 (43) 1219 (75) 1109 (35) 1082 (73) 1000 (5) 978 (8) 952 (1) 800 (9) 798 (40) 762 (8) 693 (9) 589 (1) 599 (82) 551 (3) 498 (9) 477 (6) 436 (7) 324 (14) 167 (150) 204 (13) 172 (62)

3627 (41) 3592 (73) 3492 (41) 3133 (10) 3134 (4) 1619 (538) 1601 (57) 1572 (195) 1495 (39) 1428 (298) 1367 (45) 1326 (262) 1283 (3) 1211 (31) 1103 (27) 1070 (58) 978 (10) 962 (21) 948 (1) 774 (14) 771 (37) 752 (0) 673 (4) 579 (1) 557 (86) 549 (3) 485 (12) 466 (4) 431 (11) 323 (12) 260 (229) 202 (17) 173 (1)

3645 (45) 3621 (73) 3508 (44) 3146 (9) 3120(18) 1626 (586) 1605 (39) 1579 (177) 1504 (24) 1437 (323) 1378 (45) 1334 (239) 1292 (3) 1213 (43) 1108 (21) 1072 (62) 980 (10) 961 (17) 947 (1) 779 (9) 772 (43) 755 (1) 679 (4) 580 (1) 565 (84) 546 (3) 485 (13) 466 (5) 428 (11) 323 (12) 269 (232) 200 (16) 170 (1)

3565 (105) 3591 (182) 3445 (132)

~

V I N8H9 N8H10 s ~2 0 7 H s

ooph

N8H9 N8H10 s in ph C5H s UJ C6H s Y6 C5C6 N3C2 s u7 N8H2 sciss US C4C5 N3C4 s u9 N3C2 C 2 0 s, N8H2 sciss Y I O C 2 0 N1C2 S, C6H b U I IC5C6 C4N8 s, C5H b, ring def ~ 1C 2 2 0 S, OH C6H b ~ 1 N1C6 3 C4N8 S ~ 1 OH 4 b, C4N8 s UIJ C5C6 C 2 0 S, C5H b ~ 1 OH 6 b, C5H C6H b ooph ring breathing i q 8 ring s, NH2 wag, C4C5 s ~ 1 C5H 9 C6H b ooph 00p1 u20 ring def, C 2 0 C6H s u21 ring wag oopl v22 C5H b oopl, ring wag oopl ~ 2 CO 3 C4N8 wag oopl ~ 2 OH 4 wag oopl ~ 2 ring 5 def U26 ring def u27 N8H2 tors, ring def oopl u28 N8H2 tors, ring def oopl u29 ring def oopl, N8H2 tors ~ 3 C4N8 0 C20 b u3l N8H2 b 00pl ~ 3 C4N3 2 C4N8 CO wag oopl u33 ring def oopl u3

uq

1623 (597) 1600 (63) 1569 (179) 1496 (26) 1439 (581) 1379 (74) 1320 (1 11) 1257 (16) 1193 (126) 1110 (29) 1084 (50) 980 (10) 955 (6) 782 (16) 807 (58) 751 (2) 717 (16) 601 (10) 520 (66) 569 (6) 507 (37) 451 (6) 350 (8) 297 (35)

0 b = bending, s = stretching, wag = wagging, sciss = scissoring, def = deformation, tors = torsion, in ph = in phase, ooph = out of phase, oopl = out of ring plane. Reference 43.

I5

- DFTNLZ

1

- DFT-NLZ

1

EXP

10 -

15

.....

i

I

!

1

. >.

:

,:A

.. 500

1000

1500

2000

2500

3000

3500

500

4000

1000

2). The CPU time necessary to perform 10 SCF iterations and to evaluate the gradients is given for the LDA and the gradientcorrected density functionals. Typically, we need 10-20 iterations to reach convergence in the energy (convergence criteria lo-* au) .

Acknowledgment. This investigation was supported by the “Regione Autonoma della Sardegna”. Supplementary Material Available: Structural parameters and vibrational frequencies for the U2, U3, U4, U5,and U6 tautomers of uracil and for the C2, C4, C5, and C6 tautomers of cytosine (18 pages). Ordering information is given on any current masthead page.

2500

2000

3000

3500

4000

Frequency (cm-1)

Frequency (cm-I)

Figure 4. Aminwxo tautomer of cytosine (Cl) infrared spectrum; comparison of DFT-NL2 results with experimental values.

1500

I

Figure 5. Aminc-hydroxy tautomer of cytosine (C3) infrared spectrum; comparison of DFT NL2 results with experimental values.

TABLE 1 0 Computational Expensive Associated with Uracil and Cytosine Dm Calculations (in seconds of CPU on an IBM RS6000-560 Workstation) SCF (10 iterations) gradients

U(LDA)

U(NL)

C(LDA)

C(NL)

580 110

770 215

650 130

850 230

References and Notes (1) Kwiatkowski, J. S.;Bartlett, R.J.; Person, W. B. J . Am. Chem. SOC. 1988, 110, 2353.

( 2 ) Leszcynski, J. J . Phys. Chem. 1992, 96, 1649. (3) Cieplak, P.; Bash, P.;Singh, U. C.;Kollman, P. A.J. Am. Chem. SOC. 1987, 109, 6283.

(4) Les, A.; Adamowicz, L. J . Phys. Chem. 1989, 93, 7078. (5) Les, A.; Adamowicz, L.J . Phys. Chem. 1990, 94, 7021.

5660 The Journal of Physical Chemistry, Vol. 98, No. 22, 1994

Estrin et al.

(6) Boughton, J. W.; Pulay, P. Int. J. Quantum Chem. 1993, 47, 49. (28) Perdew, J. P. Phys. Rev. 1986, B33, 8822. (7) Kwiatkowski, J. S.;Zielinski, T. J.; Rein, R. Adv. Quantum Chem. (29) Becke, A. D. Phys. Rev. 1988, A38, 3098. 1986, 18, 85. (30) Person, W. B., Zerbi, G., Eds. Vibrationallntensities inlnfraredand (8) Beak, P.; White, J. M. J. Am. Chem. SOC.1982, 104, 7073. Raman Spectroscopy; Elsevier: Amsterdam, 1982. (9) Szczesniak, M.; Nowak, M. J.;Szczepaniak, K. J . MoLStruct. 1984, (31) Scanlan, M.; Hillier, I. H. J . Am. Chem. SOC.1984, 106, 3737. 116, 387. (32) Voet, D.; Rich, A. Prog. Nucleic Acid Res. Mol. Biol. 1970,10, 196. (10) Lowdin, P. 0. Rev. Mod. Phys. 1963, 35, 724. (33) Ferenczy, G.; Harsanyi, L.; Rozsondai, B.; Hargittai, I. J . Mol.Struct. (11) Topal, M. D.; Fresco, J. R. Nature 1976, 263, 285. 1986, 140, 71. (12) Hehre, W. J.; Radom, L.; Schleyer, P. V.; Pople, J. A. Ab Initio (34) Andzelm, J. In Theory and Applications of Density Functional Molecular Orbital Theory; Wiley: New York, 1986. Approaches to Chemistry; Labanowski, J., Andzelm, J. Eds.; Springer(13) Wimmer, E.; Krakauer, H.; Weinert, M.; Freeman, A. J. Phys. Reo. Verlag: Berlin 1990. 1982, B24, 864. (35) Estrin, D.; Paglieri, L.; Corongiu, G.Unpublished results. (14) Jones, R. 0.;Gunnarsson 0. Reo. Mod. Phys. 1989, 61, 689. (36) Jasien, P. G.; Fitzgerald, G. J . Chem. Phys. 1990, 93, 2554. (15) Salahub. D. R.: Fournier. R.: Mlvnarski. P.: Paoai. I.:St. Amant. A,: (37) Brown, R. B.; Godfrey, P. D.; McNaughton, D.; Pierlot, A. P. J. Am. Ushio,'J. In Theory and Applications oiDensity Funciionai Approaches to Chem. Soc. 1988, 110, 2329. Chemistry; Labanowski, J., Andzelm, J., Eds.;Springer-Verlag: Berlin, 1990. (38) Kulakowska,I.;Geller, M.; Lesyng, B.; Wierzchowski,K.L.; Bolewska, (16) Sosa,C.;Andzelm,J.;Elkin,B.C.;Wimmer,E.;Dobbs,K.D.;Dixon, K. Biochim. Biophys. Acta 1975, 407, 420. D. A. J. Phys. Chem. 1992, 96, 6630. (39) (a) Szczesniak,M.; Nowak, M. J.; Rostkowska, H.; Szczepaniak,K.; (17) Handy, N. C.; Murray, C. W.; Amos, R. D. J. Phys. Chem. 1993, Person, W. B.; Shugar, D. J. Am. Chem. SOC.1983,105,5969. (b) Chin, S.; 97, 4392. Scott, I.; Szczepaniak, K.; Person, W. B. J. Am. Chem. SOC.1984,106,3415. (18) Berces, A.; Ziegler, T. J. Chem. Phys. 1993, 98, 4793. (40) (a) Graindourze, M.; Smets, J.; Zeegers-Huyskens, T.; Maes, G. J . (19) Estrin, D. A.; Corongiu, G.; Clementi, E. In METECC94, Methods Mol. Sfruct. 1990, 222, 345. (b) Barnes, A. J.; Stuckey, M. A. LeGall, L. and Techniquesin Computational Chemistry;Clementi, E., Ed.;Stef Cagliari, Spectrochim. Acta, Part A 1984, 40, 419. 1993; Vol. B. (41) Gould, I. R.; Vincent, M. A.; Hillier, I. H. J . Chem. Soc., Perkin (20) Kohn, W.; Sham, L. J. Phys. Reo. 1965, A140, 1133. Trans. 1992, 2, 69. (21) Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. J. Chem. Phys. 1979, (42) Vercauteren, D. P.; Paddeu, G.; Vanderveken, D. J.; Baudoux, G.; 71, 3396, 4993. Dory, M.; Clementi, E. In METECC 94, Methods and Techniques in (22) Becke, A. D. J . Chem. Phys. 1988, 88, 1053. Computational Chemistry; Clementi, E., Ed.; Stef Cagliari, 1993; Vol. B. (23) Sim, F.; Salahub, D. R.; Chin, S.; Dupuis, M. J . Chem. Phys. 1991, (43) Szczesniak, M.; Szczepaniak, K.; Kwiatkowski, J. S.;KuBulat, K.; 95, 43 17. Person, W. B. J. Am. Chem. Soc. 1988,110, 8319. (24) Sim, F.; %.-Amant, A.; Papai, 1.; Salahub, D. R. J. Am. Chem. SOC. (44) Les, A.; Adamowicz, L.; Bartlett, R. J. J. Phys. Chem. 1989, 93, 1992, 114, 4391. 400 1. (25) Gunnarsson, 0.;Lundqvist, B. I. Phys. Reo. 1976, B13, 4274. (45) Leszcynski, J. Int. J . Quantum Chem., Quantum Biol. Symp. 1992, (26) Vosko, S.H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200. 19, 43. (27) Perdew, J. P. Phys. Reo. 1986,833,8800; 1986,834,7406(erratum). (46) Gould, I. R.; Hillier, I. H. Chem. Phys. Lett. 1989, 161, 185.