J. Phys. Chem. 1995, 99, 11417- 11423
11417
Effect of the Reference Geometry and the Exchange Correlation Functional on the Vibrational Frequencies Calculated by Density Functional Methods. The Examples of Benzene and Nickel, Chromium, and Iron Carbonyls Attila Berces and Tom Ziegler" Department of Chemistry, The University of Calgary, Calgary, Alberta E N 1N4, Canada Received: March 24, 1995; In Final Form: May 15, 1995@
The vibrational frequencies of benzene, Ni(C0)4, Cr(C0)6, and Fe(C0)S have been calculated with the local density approximation (LDA) at various reference geometries. The reference geometries were taken from experiment and from geometry optimizations based on the LDA method as well as LDA augmented with nonlocal exchange-correlation corrections (LDA/NL). We found that the LDA method reproduces the observed harmonic frequencies remarkably well at the experimental reference geometry and most of the error at the LDA optimized structure can be related to reference geometry effects. For benzene, Cr(C0)6, and Ni(C0)4, the average deviations (and average percentage deviations) were 15.3 cm-' (1.32%), 5.1 cm-' (1.55%), and 7.0 cm-I (2.34%) at the experimental reference geometry, respectively. At the LDA optimized geometry, the average deviations were 32.6 cm-' (1.80%), 29.4 cm-' (7.56%), and 33.9 cm-' (10.6%), respectively. For Fe(C0)5, the best agreement was found at the LDA/NL optimized reference geometry rather than the experimental structure, and this was attributed to significant uncertainties in the experimental structure. A review is provided of the frequency assignments for the iron pentacarbonyl spectra based on the calculated results, and some alternative assignments are suggested. We found no significant difference between the vibrational frequencies of the metal carbonyls calculated by local and nonlocal DFT methods, provided the reference geometries were the same.
I. Introduction A number of recent studies have evaluated density functional theory (DFT) as a practical tool for frequency calculations, including our own investigations in which special emphasis was put on transition metal complexes. The general conclusion seems to be that the current approximate density functionals afford quite adequate estimates of frequencies, even at the simplest level represented by the local density approximation (LDA). However, the sources of deviation between calculated and observed frequencies are not well understood. Most studies in this field compare frequencies obtained from calculations with different exchange correlation functionals.I It is generally observed that the gradient corrected (or nonlocal) density functional methods provide more accurate geometries and vibrational frequencies than the simple local DFT methods. The improved geometry at the higher level of calculation is especially remarkable for the transition metal complexes. The more accurate frequencies at the nonlocal level can be explained partially by the improved reference geometry and partly by the contribution of the nonlocal gradient corrections to the energy Hessian. In the present study we look at these two effects separately by performing frequency calculations by LDA method at various reference geometries. Further, we compare frequencies obtained at the same reference geometry by the two different levels of DFT calculations. In quantum mechanical calculations of vibrational frequencies generally the optimized geometry is utilized as a reference geometry. This choice greatly simplifies the calculations since the frequencies can be obtained by a simple diagonalization of the mass-weighted Cartesian force constants. Most major quantum mechanical programs do not offer other choices of reference geometry. Any geometry other than the optimized 'Abstract published in Advance ACS Absrructs, July 1, 1995.
0022-365419512099-11417$09.0010
geometry would introduce nonzero forces on the atoms at the reference point, and methods that circumvent the nonzero force dilemma are not generally implemented. However, we shall in our frequency calculations at structures other than the optimized geometry take the nonzero forces into account in a proper way. The use of observed structures as a reference was pioneered by Schwendeman,2who observed that the calculated frequencies of diatomic molecules improved when they were evaluated at the experimental geometry. For polyatomic molecules, Blom and Altona introduced empirical correction tenns to account for the deficiencies of the Hartree-Fock optimized geometry and suggested the empirically corrected geometry as a refere n ~ e .Pulay ~ and Fogarasi also pointed out the importance of the reference geometry and further developed the original correction method due to Blom and Altona by introducing offset forces as correction factor^.^ Pulay et al. also pointed out that a change of only 0.016 8, in the bond length produces 10% error in the stretching force constants due to anharmonicities, based on a simple estimate from a Morse p~tential.~More recently, Allen and CsBsz6r investigated the effect of the reference geometry on the harmonic as well as anharmonic force constants.6 We have briefly discussed how different geometries influenced the frequencies of benzene calculated by DFT based methods in a previous study.' However, the nonzero forces were not treated properly in that investigation. The question of reference geometry otherwise has not been studied in connection with density functional methods. The present study will cover the vibrational frequencies of benzene as well as the metal carbonyls Cr(C0)6, Fe(C0)5, and Ni(C0)4 calculated by density functional methods at several reference geometries. The reference geometries were obtained from geometry optimization using local (LDA) and nonlocal exchange correlation function0 1995 American Chemical Society
11418 J. Phys. Chem., Vol. 99, No. 29, 1995 als, and from experiment. We also compare the performance of the LDA and LDA/NL methods at a given reference geometry.
II. Computational Details The reported calculations employed the ADF program system developed by Baerends et aLgaand vectorized by Ravenek.8bte Velde et al? developed the numerical integration procedure. The atomic orbitals on the metal centers were described by an uncontracted triple-5 STO basis set,I0 while a double-g STO basis set was used for C, 0, and H. The basis sets on H, C, = l.O), 3d = 2.5), and and 0 were augmented with 2p = 2.5) polarization functions, respectively. The ls2 3d configuration on carbon and oxygen as well as the ls22s22p6 configuration of the metals was assigned to the core and treated by the frozen-core approximation.8a A set of auxiliary” s, p, d, f, and g STO functions, centered on all nuclei, was used to fit the molecular density and represent the Coulomb and exchange potentials accurately in each SCF cycle. Optimized geometries in this study were calculated based on the energy expression from the local density approximation’2a (LDA) augmented by nonlocal corrections to exchange12band correlation.12c We shall refer to this method as the LDA/NL scheme, which we have implemented self-c~nsistently~~ in energy and gradientI4calculations. The geometry was optimized based on the GDIIS technique.I5 The Cartesian force constants and dipole moment derivatives were calculated by numerical differentiation of the energy gradients16 using Cartesian displacements. We have implemented an automatic scheme for the transformation of symmetry related Cartesian force constants.Ik Therefore, only the symmetry unique displacements had to be made. We have used a sufficient number of integration points NPto ensure a numerical accuracy of 1.0 cm-I for most vibrational frequencies. The numerical accuracy was tested by observing the stability of the calculated frequencies as a function of NP. For frequencies evaluated at other points than the optimized geometry the forces or energy gradients are different from zero. Therefore, when the force constants are transformed into internal coordinates, q, terms due to the gradient also have to be taken into account.6 The force constants in internal coordinates can be expressed in terms of the Cartesian energy gradients and energy Hessians as
(cyd
(cyp
(c:,
We have implemented this transformation into a locally developed program packageI7 based on Schachtschneider’sforce field program.I8 This program was used for the normal coordinate analysis. The frequencies obtained by this technique are not completely independent of the selection of intemal coordinate^.^^,^ In practice, a physically meaningful set, like the natural intemal coordinates suggested by Fogarasi et al. works well.& For a detailed discussion on the meaning of the correction term the reader is referred to a recent study by Zhou et a1.4d
III. Experimental Reference Geometries In force field and frequency calculations, ideally the equilibrium (re)geometries should be used as a reference point. Care should be taken when experimental geometries are used, since there might be significant uncertainties associated with the experimental values due to the approximations in the models
BCrces and Ziegler used in the evaluation. Empirically, the equilibrium structure is available only for diatomic molecules and for the simplest polyatomic molecules. The geometrical parameters obtained by experiment correspond to vibrating and rotating molecules, which can also interact with other molecules in condensed phases. Accordingly, the experimental structures determined under differennt conditions and by different experimental methods carry different information. Further, as more accurate geometries are desired, the representation of empirical geometry becomes increasingly important. The rg structures is the thermal average of interatomic distances determined by diffraction methods, which is in most cases close to the equilibrium structure. From this respect it is an advantage that empirical rg structures are available for all transition metal carbonyls studied here. Jost and co-workers determined the rg geometry of chromium hexacarbonyl from a neutron diffraction study in the solid state.I9 The solid state geometry can be substantially different from the gas phase geometry due to intermolecular interactions. These interactions, however, mainly affect the bond angles and torsion angles, especially the soft degrees of freedom. The average bond length should not be significantly different in the gas and solid states. Jost et al. determined the rg structural parameters based on two models for the correction of vibrational averaging. One of the two models (expz in Table 2) resulted in better agreement with a previous X-ray study by Whitaker and JefferySz0We included both sets of parameters in our study, since we suspected that the experimentalists had rejected the correct geometry (expl in Table 1). We retum to this point in the discussion. The experimental geometry for Ni(CO)4 was determined from room temperature gas phase electron diffraction experiments.21 The observed values are corrected for vibrational motions, and the corrected rg structure is used as reference geometry in our study. This experimental geometry is the most ideal for our calculations, since it is expected to yield results very close to the gas phase equilibrium structure. Most recently Beagley and Schmidling derived rg structural parameters of iron pentacarbonyl from electron diffraction experiments utilizing an approximate harmonic force field to correct for vibrational effects.22 The same study reported rg parameters corresponding to two previous experimental meas u r e m e n t ~ . The ~ ~ averages of the axial and equatorial bond distances, and the average difference between the equatorial and axial Fe-C bonds, were determined. The experimental structures suggest that the equatorial FeC bond is slightly longer than the axial FeC bond, which conflicts with the theoretically predicted structures at both the LDA and the LDA/NL level. Further, the axial and equatorial CO bonds are assumed to be equal in the experimental work, which is a poor approximation for our study. In the two earlier studies23the average CO bond length was found substantially shorter, 1.146 and 1.148 A, than in the most recent one. The difference between the axial and equatorial FeC bond length varies between 0.012 and 0.026 A, depending on which experimental structure we take. Due to these uncertainties, these experimental geometries are not an appropriate reference geometry for frequency calculations. The most recent experimental geometry is included in our study along with the optimized structures for completeness. For benzene, there is no experimental average structure available, the ro and r, structures are determined. The ro (expl in Table 1) structure was determined by S t ~ i c h e f f . The ~~~ advantage of the ro structures is that the rotational constants can be measured very accurately, and there is no uncertainty due to different assumptions. The r,,, “mass dependence” structure represents approximate equilibrium structure based on an approximate relationship
Vibrational Frequencies of Benzene and Metal Carbonyls
TABLE 1: Vibrational Frequencies (cm-') of Benzene Calculated at Different Reference Geometries ref geom
CWA
A2, Bi, Bzu
Ezg
El,
Azu B2g El,
Ezu
2 1 3 13 12 14 15 7 8 9 6 20 19 18 11 5 4 10 17 16
LDA calculations expZb LDA
PFBc
10.84 1.397
1.0862 1.3902
1.094 1.388
1.077 1.395
expd 3191 994.4 1367 3174 1010 1309.4 1149.7 3174 1607 1177.8 607.8 3181.1 1494 1038.3 674.0 990 707 847.1 967 398
3193.6 982.7 1321.5 3159.8 983.5 1336.4 1126.4 3170.2 1579.5 1150.3 594.0 3184.4 1453.1 1029.5 662.9 990.9 706.5 834.2 959.8 400.3
3173.6 995.9 1322.6 3139.0 989.7 1362.2 1129.3 3149.6 1598.7 1154.0 598.9 3164.1 1461.7 1037.5 661.8 992.6 712.0 835.4 962.9 404.2
3101 1004 1314 3065 993 1379 1125 3075 1610 1150 602 3091 1462 1039 664 985 713 830 952 399
3257.8 984.5 1329.9 3224.5 984.8 1350.8 1145.6 3237.5 1591.4 1170.1 589.4 3249.0 1463.5 1034.7 689.8 1038.5 715.2 819.3 998.5 434.8
av dev av % dev
15.31 1.32
16.22 1.25
32.56 1.80
30.47 2.52
CC/A
AI,
expi"
Reference 24a. Determined using the rotational constants of four molecules (ref 24b) Approximation based on empically corrected ab initio values. Reference 26. Reference 29.
D6h
between the rotational constants for the ro structure of the parent molecule and its isotopic substitution^.^^ This method is usually not very accurate for determination of the positions of light atoms, especially hydrogen. The r, structure of benzene (exp2 in Table 1) was reported by Pliva et ~ 2 1 . ~ ~ ~ Pulay et al. derived theoretical approximations to the equilibrium bond distances of benzene from an empirically corrected Hartree-Fock calculation.26 This estimated structure (PFB in Table 1) is also used in our study as a reference geometry. Pulay et al. pointed out in their more recent publication that this CH reference bond length for benzene might be somewhat too short.*' The equilibrium CC and CH bond length is expected to be shorter than the ro values and the correction should be larger for the CH bond. Indeed, the PFB geometry is consistent with this expectation. On the other hand, the r, CH bond is longer than the ro value, which also points out the weakness of the r, method for light atoms. Further, the vibrational averaging for the CC bond (Le., ro-r,,,) seems somewhat too large. Tentatively, the CC bond of the PFB geometry seems the best approximation for the equilibrium at 1.395 A. Estimates based on theoretical arguments and experimental data by Handy et al. put the equilibrium CH bond length between 1.0800 and 1.0825 A, which also seems reasonable.28
IV. Comparison of Vibrational Frequencies A. Comparison of LDA Frequencies at Different Reference Geometries. In this section we study the differences in the vibrational frequencies calculated at various reference geometries. The calculated frequencies of benzene and chromium, nickel, and iron carbonyls at experimental, LDA and LDA/NL optimized reference geometries are compared to the observed vibrational frequencies in Tables 1-4. All the theoretical data refer to simple LDA force field calculations. Benzene. The experimental vibrational f r e q ~ e n c i e sof ~~ benzene listed in Table 1 represent harmonic modes, although the anharmonic corrections of certain vibrations have been
J. Phys. Chem., Vol. 99, No. 29, 1995 11419 questions by Handy et aLZ8 The CH stretching modes of benzene, vn,with n = 2,7, 13, and 20, are very good candidates for studying the effect of the reference geometry on calculated frequencies. These vibrations are highly sensitive to the applied CH bond length, Table 1, and they are not coupled to other vibrations. While the CH stretching frequencies at the LDA reference geometry are underestimated by about 100 cm-', they are systematically above the experimental values by 60 cm-l at the PFB geometry. The overestimation of the CH frequencies at the PFB geometry is likely related to the uncertainty associated with the empirical correction term used to obtain that CH bond length. (See previous section.) The first experimental geometry, expl of Table 1, yields the best results for the CH stretching frequencies, with an average deviation of only 5.2 cm-I. Although this agreement between experiment and theory is exceptionally good, one has to be careful with the interpretation. This reference geometry does not reflect the equilibrium value. Therefore, the excellent agreement between theory and experiment is related to cancellation of errors to a certain extent. On the basis of the approximate CH bond length by Handy et al., which is between 1.0800 and 1.0825 A, one can estimate what the LDA CH stretching frequencies would be at the equilibrium geometry. These estimates show that the LDA CH stretching frequencies would be at most 30 cm-' different from the experimental values, on average. Nonetheless, most of the 100 cm-' error at the LDA reference geometry can be accounted for by the too long CH bond length. The CC stretching vibrations, v,, with n = 1, 8, 14, and 19, are also sensitive to the reference geometry; however, for benzene these are not as straightforward to interpret as the CH stretching vibrations. One characteristic of the benzene force field is the strong coupling between the CC bond stretchings. For certain symmetries the CC coupling constants have primary importance, while couplings with CH rock and ring deformation coordinates lead to strongly mixed vibrations. Note that the empirical anharmonic corrections of these frequencies are not as well established as for the CH vibrations, which might also contribute to the increased relative error, compared to the CH stretching frequencies. These values, however, also indicate that some of the interaction force constants are not calculated accurately by LDA method, since the calculated CC stretching frequencies show nonsystematic error at all reference geometries. Bending frequencies are less sensitive to the change in the bond length. They change mainly because they are mixed with the CC stretching vibrations. Cr(COk. The calculated vibrational frequencies of Cr(C0)6 are listed and compared to observed frequencies in Table 2. The empirical frequencies were determined by Jones et aL30 The experimental CO stretching frequencies are corrected for anharmonicity, while other frequencies are the observed fundamentals. The neglected corrections are expected to be less than a few cm-I. The CO stretching frequencies of metal carbonyls are sensitive not only to the CO distance but also to the metalcarbon distance due to strong coupling. The geometry labeled expl yields the CO stretching frequencies in excellent agreement with the experiment (9.6 cm-I or 0.5% deviation). The LDA predicted and experimental (expl) CO bond lengths are the same; therefore, the calculated stretching frequencies are also similar at these two reference geometries. However, the seriously underestimated CrC bond length by the LDA method slightly increases the CO stretching frequencies through coupling with CrC stretching vibrations. The LDA/NL method yields more accurate CrC bond distances, while it overestimates the CO
BCrces and Ziegler
11420 J. Phys. Chem., Vol. 99, No. 29, I995 TABLE 2: Vibrational Frequencies (cm-') of Chromium Hexacarbonyl Calculated at Different Reference Geometries ref LDA calculations geom explb e ~ p 2 ~LDA/NL LDA 1.918 1.917 1.862 CCIA 1.916 COIA
1.147
1.141
1.154
1.147
2139.2 379.2 2045.2 390.6 364.1 2043.7 668.1 440.5 97.2 532.1 89.7 510.9 67.9
2141.4 383.3 2050.1 390.6 359.3 203 1.5 687.2 441.3 100.6 535.3 90.7 513.1 63.9
2179.0 385.2 2090.9 393.9 364.9 207 1.O 690.9 448.4 95.7 538.8 90.6 516.6 67.7
2093.1 381.2 2001.0 388.4 356.0 1982.0 683.7 438.2 99.5 531.9 89.9 509.5 63.3
2159.4 440.1 2065.2 447.1 380.0 2044.9 747.2 504.4 111.3 568.3 98.6 536.8 67.8
av dev av % dev
5.10 1.55
11.07 1.38
13.03 1.86
29.39 7.56
expa Ai,
E Fig
FI,
F2, Fzu
Reference 30. Reference 19. distance by 0.007 A. The longer CO bond results in 50 cm-' downshift for the calculated CO stretching frequencies. The second experimental geometry (expz) yields the CO stretching frequencies with substantially larger error (35.5 cm-', 1.7%). An analysis of this error is provided in the discussion. While the LDA/NL reference geometry yields the worst CO stretching frequencies, the CrC frequencies are calculated most accurately at this reference geometry. There is, however, no significant difference between the predictions of CrC stretching frequencies at exp 1, exp2, and LDA/NL reference geometries; the average deviations are 7.5 cm-' (1.8%), 6.1 cm-' (1.4%), and 5.0 cm-' (1.2%), respectively. The LDA method underestimates the CrC bond distance; as a result the corresponding CrC stretching frequencies are overestimated by as much as 64 cm-' (15%). The deformation frequencies are similar at the two experimental and the LDA/NL reference geometries, and they are all in excellent agreement with the experimental values. Although stretching frequencies are generally more sensitive to the change in reference bond length, large deviations in the bond lengths affect the bending frequencies as well. This is apparent from the calculated values at the LDA reference geometry. Ni(C0)d. The calculated and experimental vibrational frequencies for Ni(C0)d are listed in Table 3. The empirical frequencies in Table 3 were determined by Hedberg et aL2' The CO stretching modes represent harmonic frequencies, while the rest of the frequencies are not corrected for anharmonicity. The calculated CO stretching frequencies at the experimental geometry are within 1 and 7 cm-', compared to the observed harmonic values. Since the LDA value of the CO bond length is close to the experimental, the corresponding calculated frequencies are also in good agreement with experiment. The LDA/NL CO bond length is 0.009 A too long which results in 60-70 cm-' underestimation of the CO stretching frequencies. While the LDA reference geometry is appropriate for the CO stretching frequency, it seriously overestimates the NiC stretching frequencies due to the underestimated NiC bond length. These observations for both the CO and metal-carbon stretching frequencies are in line with those observed for Cr(C0)6. There is one experimental vibrational frequency that is incompatible with the present calculated data. The 380 cm-' frequency in the E symmetry block is most likely incorrectly determined, as previously pointed out by Sosa et aL3'
TABLE 3: Vibrational Frequencies (cm-') of Nickel Tetracarbonyl Calculated at Diaerent Reference Geometries ref LDA calculations geom exp" LDA/NL LDA NiCIA 1.838 1.844 1.781 COIA 1.141 1.150 1.143 expa A
E FI
K
2154.1 370.8 380.W 62.0 300.0 2092.2 458.9 423.1 79.0b
2153.0 369.1 465.6 58.7 278.9 2084.9 459.8 428.1 67.2
2086.3 361.8 442.8 38.4 257.2 2017.2 45 1.6 419.7 56.3
2157.6 433.4 485.3 47.0 283.9 2088.5 518.9 491.9 67.3
av dev av 8 dev
7.02 2.34
3 1.84 8.90
33.87 10.63
Reference 21. Not included in the calculation of average deviations, see text.
TABLE 4: Vibrational Frequencies (cm-I) of Iron Pentacarbonyl Calculated at Different Reference Geometries LDA calculations ref geom expb LDA/NL LDA act. axial FeCIA 1.810 1.819 1.769 COIA 1.153 1.152 1.145 equatorial FeC/A 1.831 1.817 1.766 CO/A 1.153 1.156 1.149 expa Al'
A*' E'
AT
E"
2143 2064 443 413 383' 2038 645 488 429 104 64' 2064 619 474 100 543 375 97
2101.4 2020.2 453.3 407.7 351.9 2007.3 658.0 466.6 425.8 98.9 49.6 2010.0 616.7 486.4 101.3 550.8 365.9 93.5
2088.6 2006.7 441.3 416.5 352.6 1984.4 656.5 476.2 428.8 97.9 47.0 2007.5 613.0 476.8 102.1 550.2 365.6 92.6
2158.0 2074.8 500.3 476.4 373.1 2052.1 713.2 533.4 453.9 104.3 49.4 2072.9 680.3 530.3 108.2 581.3 387.5 98.1
av % dev
14.94 2.36
16.36 2.30
28.78 5.68
R
ia IR, R
IR
R
aReferences 32 and 33. bReference 22. 'Not included in the calculation of average deviations, due to uncertainty in the assignments. Fe(C0)j. In Table 4, we compared the observed and calculated vibrational frequencies of iron carbonyl at three different reference geometries. The three reference geometries are the optimized LDA, LDNNL, and the experimental structures. The observed CO stretching frequencies listed in Table 4 include anharmonic corrections, while the rest of the frequencies are the observed fundamentals. The assignments in Table 4 are based on several experimetnal results. Before we compare frequencies at different reference geometries, we discuss the frequency assignments. Empirical frequency assignments for iron pentacarbonyl differ significantly by different authors. We found that our calculations c o n f i i e d most of the assignments by Jones et al.,32while more recent assignments are inconsistent with the calculated results.33 Most studies, including the one by Jones et al., agree that the observed absorption frequency at 542 cm-' is of E' symmetry, since it is observed in both the IR and the Raman
J. Phys. Chem., Vol. 99,No. 29, 1995 11421
Vibrational Frequencies of Benzene and Metal Carbonyls spectra. However, this assignment is not consistent with our calculations. In light of the good agreement between calculations and observations for similar frequencies of nickel and chromium carbonyls, we suggest that the observed IR band is probably not a fundamental, and the Raman band is therefore not E’ symmetry. The low intensity of this observed IR band, and inconsistencies in the isotopic shifts32 of this frequency, further support this assignment. Therefore, only the observed Raman line is a fundamental, which allows one to assign this frequency to the E” symmetry. This assignment is consistent with the calculated value of 550 cm-’. Two frequencies around 470-480 cm-I of E’ and A2” symmetries, respectively, are expected based on the calculations. We found these in the IR spectrum; however, one of these was assigned to be a combination band. We modified the assignment of the observed 474 cm-’ IR band from E’ to A*’’ symmetry, which allows the 488 cm-’ observed Raman band (originally E”) to be assigned to the E’ symmetry. Further, the observed 429 cm-’ IR line can be assigned to A; symmetry, in line with the observation that this band does not have a counterpart in the Raman spectrum. The proposed assignments are consistent with the selection rules and they are in good agreement with the theoretical predictions. The CO stretching frequencies of both chromium and nickel carbonyls were calculated with a few cm-I accuracy at the empirical reference geometry. For Fe(C0)s however, such accuracy cannot be obtained due to ambiguities related to the empirical CO bond length. The experimental geometry was determined assuming equal axial and equatorial CO bond lengths, whereas the calculation predicts them to be different by 0.004 A. This 0.004 A difference in the CO bond length translates into about 30 cm-’ shift in the calculated CO stretching frequencies. Further, the earlier experiments predicted 0.007 8, shorter average CO bond length than the most recent one. For these reasons, we could not calculate CO stretching frequencies with high accuracy at an empirical geometry. Except for the results at empirical reference geometry, the Fe(C0)5 results are similar to those for chromium and nickel carbonyls. The metal-carbon frequencies significantly improved by taking the LDAlNL optimized or the experimental geometry as reference. The CO stretching frequencies of iron pentacarbonyl agree the best with the results at the LDA reference geometry. Deformation frequencies were similar at all geometries. The only notable exceptions are the MCO bendings in the 400-500 cm-I region at LDA geometry, which were increased by the too short MC bond length. B. Comparison of LDA and LDA/NL Vibrational Frequencies at a Given Reference Geometry. In this section we compare the vibrational frequencies of Cr(C0)6, and Fe(C0)s calculated by LDA and LDAlNL methods at experimental reference geometries. The frequencies reported here are always tested for numerical stability with respect to the number of integration points, and about 1 cm-’ numerical accuracy is ensured for most frequencies. We have experienced difficulties with the Ni(C0)4 LDA/NL calculation, where reasonable accuracy could not be reached for some low-frequency vibrations. The molecular orientation in the Cartesian coordinate system plays a role in the signal/noise ratio when force constants are evaluated by Cartesian displacements. For Cr(C0)6, for example, it is possible to orient the molecule in such a fashion that all Cartesian displacements are either perpendicular or parallel to the bonds. Such orientation is not possible for the tetrahedral Ni(C0)4, which explains the unique numerical difficulties. We have taken the experimental geometry as a reference geometry for the comparison of LDA and LDA/NL frequencies.
TABLE 5: Com arison of LDA and LDA/NL Vibrational Freauencies (cm- ) at Experimental Geometrv
P
AI, E
FI, Flu
Fig Fzu
2141.4 383.3 2050.1 390.6 359.3 2031.5 687.2 441.3 100.6 535.3 90.7 513.1 63.9
2140.6 389.3 2046.3 393.8 368.4 2028.1 697.9 451.6 108.4 548.3 92.8 525.8 64.8
AI’
All E’
AT
E“
2101.4 2020.2 453.3 407.7 351.9 2007.3 658.0 466.6 425.8 98.9 49.6 2010.0 616.7 486.4 101.3 550.8 365.9 93.5
2104.9 2019.4 455.4 410.9 362.5 2008.4 668.2 474.2 435.9 106.4 44.1 2010.7 627.4 493.3 110.8 565.6 376.0 95.6
TABLE 6: Comparison of Vibrational Frequencies of Cr(C0)a Obtained by Truncated and Complete Force Field Transformation ref geom (LDA/NL) complete truncated AI, 2093.1 2093.1 381.2 381.2 E 2001.0 2001.0 388.4 388.4 FI, 356.0 409.1 Flu 1982.0 1982.0 683.7 710.5 438.2 441.8 99.5 139.9 F2, 531.9 567.0 89.9 119.3 F2u 509.5 547.1 63.3 95.6 This comparison could have been done at any reference geometry; the results would be similar. Table 5 includes the calculated frequencies of Cr(C0)6, and Fe(CO)S by the two levels of density functional methods. We have not included benzene in this study, since the nonlocal correction is expected to have little influence on organic molecules. While we have seen significant differences between the geometries at the LDA and LDA/NL level, for frequencies such large deviations are not seen in Table 5. The largest deviation in the stretching frequencies is the 10 cm-’ increase of the FI, CrC stretching frequency of Cr(C0)6 (687.2 cm-’ by LDA). It is especially remarkable that the MC stretching frequencies are not changed significantly by the nonlocal corrections, even though these corrections had an important effect on the MC bond length. The E” FeCO bending frequency of Fe(C0)5 increased from 551 to 566 cm-I, which represents the largest deviation for deformational frequencies. Comparing these numbers to experiment is difficult, since the deviations between LDA and LDAML results are similar magnitude to the neglected anhannonic corrections.
V. Discussion We start our discussion with some comments on the differences between the transformation of force fields at stationary and nonstationary geometries. In Table 6 we show the effect of neglecting the second term at the right hand side of eq 1 in the force field transforamtion at nonstationary points for Cr-
11422 J. Phys. Chem., Vol. 99, No. 29, 1995
(CO)6. This comparison shows that the pure stretching vibrations (e.g., all AI, and E) are not affected by the truncated transformation. For stretching force constants the second term of the transformation vanish, since the second-order derivatives of stretching coordinates with respect to Cartesians are zero. For deformation coordinates the truncated expression yields erroneous results. The error in the frequencies by the truncated transformation is sometimes half the magnitude of the frequencies themselves in the low-frequency region. When the frequencies are calculated by the diagonalization of the massweighted Cartesian force constants, the resulting frequencies must be the same as the ones obtained from the truncated transformation by definition. Consequently, mass-weighted Cartesian coordinates are not appropriate for the calculation of vibrational frequencies at geometries other than the optimized geometry. We continue our discussion with some generally observed trends for the CO and metal-carbon bond length and corresponding vibrational frequencies by DFT methods. For all three transition metal complexes, calculations at the LDA reference geometry produced CO stretching frequencies in excellent agreement with experiment. Some previous results on other organic molecules and organometallic complexes also support this general trend.34 For nickel carbonyl the LDA CO bond length differs only by 0.002 8, from the experimental rg value. For iron carbonyl the most recent experimental CO bond length is uncertain. However, our LDA calculated average CO distance of 1.147 8, is between the 1.148 and 1.146 8, values suggested by the two previous experimental studies. For chromium carbonyl, the CO distance of expl structure is in excellent agreement with the LDA results. Looking at the nonlocal CO bond lengths, they are 0.007 8, longer than the corresponding LDA values, without exception. Further, calculations by Fan and Ziegler34bsuggest that the CO bond distance of Mo(CO)6 also follows this general trend; the LDA CO bond length is in excellent agreement with the experimental values, while the LDA/NL bond length is 0.006 8, longer. The CO stretching frequencies being very close to experimental values at the LDA reference geometries also support that the error in the CO bond length is systematic, and the LDA values are probably close to the equilibrium. This trend helps us decide between the alternative experimental geometries for chromium and iron carbonyls. On the basis of these results, we believe that the geometry labeled expl for chromium carbonyl, and the earlier reported average CO bond length of Fe(CO)5 by Beagley et al.23aor Almenningen et aZ.,23bare closer to the equilibrium values than the alternative geometries. The experimental geometry exp2 for Cr(C0)6 agrees better with earlier X-ray results. While the X-ray is scattered on the electronic cloud, the neutron diffraction (ND) methods is a nuclear phenomena. Therefore, these two methods can produce slightly different results. A different general trend is observed for the metal-carbon stretching frequencies of the transition metal complexes studied here. The MC bonds are always underestimated by 0.05-0.06 A, compared to the experimental or LDA/NL values. As a consequence, the metal-carbon stretching frequencies are overestimated by 60-70 cm-I. The MC bond length calculated at the LDA/NL level always agree within a few milliangstroms with the experimental rg values. As the reference MC bond length improves, the agreement between the LDA calculated frequencies improves so that most theoretical frequencies are within a few cm-' compared to the experimental values. For the CH stretching frequencies of benzene, our observation was similar: most of the 100 cm-' error at the LDA reference
BBrces and Ziegler geometry was recovered by simply improving the reference bond length. From the trends of the CH, MC, and CO stretching frequencies at different reference geometries, one can conclude that most of the error in the stretching frequencies calculated at the LDA level at LDA-optimized geometry is related to the erroneous reference bond distance. We recall that the LDA and LDA/NL stretching frequencies were similar for the metal carbonyls, provided the reference geometry was the same. Consequently, the gradient-corrected exchange correlation functional affords significant improvement in the geometry compared to the local methods but has little effect on the calculated force constants. The previous observation that LDA/NL methods provide better frequencies than LDA methods is mainly due to the improved reference geometry. Although the reference geometry in most cases improves using the nonlocal functional, the CO bond distance becomes too long, and the corresponding stretching frequencies are too low. We also found that deformational frequencies are less sensitive to small changes in the reference bond length. On the other hand, bending frequencies sometimes differ between LDA and LDA/ NL calculations at a given reference point. These deviations, however, are still not significant. The observation that the nonlocal functional does not influence the force constants compared to LDA at a given reference point is not surprising considering experience from ab initio theory. When in ab initio methods one compares the HartreeFock method with correlated methods, the effect of the correlation on the bonding energy is very significant. The contribution of the correlation energy to various orders of energy derivatives decreases with increasing order. The correlation energy term can be as large in magnitude as the HF bonding energy itself, while the correlated geometries (Le., gradientrelated properties) differ only by a few percent and force constants by even less. M e n and C s h z k found that the quartic force fields of F20 and N20 calculated by high-level ab initio coupled cluster methods are reproduced remarkably well by the simple restricted Hartree-Fock predictions utilizing the same reference geometry.6 The electron correlation did not contribute significantly to the values of anharmonic force constants. The trend for the contribution of gradient-corrected exchange correlation functional to bonding energies, geometries, and harmonic and anharmonic force constants should be similar. Further, a detailed analysis of the effect of the reference geometry on various orders of force constants by Allen and CsBszAr showed that the major part of the error introduced by the erroneous reference geometry is related to the nuclearnuclear repulsion term.6 This effect can be shown by simply examining how the electronic and nuclear part of the total molecular energy contribute to various orders of force constants. Let us separate the electronic and nuclear potential energies as follows:
By differentiating eq 2 with respect to one of the intemal stretching coordinates RABtwice, we get
In eq 3, some terms related to the change in the interatomic distances other than RAB are not written out and they play a secondary role. It is apparent from eq 3 that this stretching force constant is sensitive to the value of RABdue to the I I R A B ~
Vibrational Frequencies of Benzene and Metal Carbonyls dependence of the second term. The nuclear repulsion term is treated by classical physics in the expression of the total molecular energy. Therefore, when one performs a calculation at a reference geometry significantly different from the (empirical) equilibrium geometry, one introduces an error in the inexpensive classical mechanical term that may ruin the accuracy of the expensive quantum mechanical calculation.
VI. Conclusion This study showed that the simple local density functional method is very accurate for the calculation of harmonic force fields and vibrational frequencies. However, accurate calculations require reference geometries that are close to the equilibrium geometry and the calculated frequencies are sensitive to the reference bond length. When LDA frequencies are calculated at the optimized geometries, most of the error in the frequencies can be clearly linked to the inaccurate reference bond length. Only the error in the highly coupled CC stretching frequencies of benzene could not be explained by reference geometry effect. Deformational force constants and frequencies are less sensitive to geometry effect. Higher level of density functional theory may improve the geometries significantly, but they do not influence the force constants compared to LDA calculation at a given reference point. Consequently, improvements of density functional frequency calculations should focus on obtaining better geometries. When different levels of theory are compared, the effect of reference geometry should be treated separately. Acknowledgment. This investigation was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) as well as the donors of the Petroleum Research Fund, administered by the American Chemical Society (ACS-PRF 27023-AC23). The Academic Computing Services of the University of Calgary is acknowledged for access to the IBM6000/RISC facilities. References and Notes (1) Delley, B.; Wnnn, M.; Luthi, H. P. J. Chem. Phys. 1994,100,5785. (2) (a) Schwendeman, R. H. J . Chem. Phys. 1966, 44, 556. (b) Schwendeman, R. H. J. Chem. Phys. 1966, 44, 2115. (3) Blom, C. E.; Altona, C. Mol. Phys. 1976, 31, 1377. (4) (a) Pulay, P.; Fogarasi, G.; Pang, F.; Boggs, J. E. J . Am. Chem. Soc. 1979, 101, 2550. (b) Forgarasi, G.; Pulay, P. In Vibrational Spectra and Structure; Dung, J. R., Ed.; Elsevier: New York, 1986; Vol. 14, pp 125-219. (c) Fogarasi, G.; Zhou, X.; Taylor, P. W.; Pulay, P. J . Am. Chem. Soc. 1992, 114, 8191. (d) Zhou, X. et al. Spectrochim. Acta 1993, 49A, 1499. ( 5 ) Pulay, P.; Lee, J. G.; Boggs, J. E. J. Chem. Phys. 1983, 79, 3382.
J. Phys. Chem., Vol. 99, No. 29, 1995 11423 (6) Allen, W. D.; CsBszL, A. J. Chem. Phys. 1993, 98, 2983. (7) Berces, A.; Ziegler, T. J. Chem. Phys. 1993, 98, 4793. (8) (a) Baerends, E. J.; Ellis, D. E.; Ros, P. Chem. Phys. 1973, 2, 41. (b) Ravenek, W. In Algorithms and Applications on Vector and Parallel Computers; te Riele, H. J. J., Dekker, Th, J., van de Vorst, H. A,, Eds.; Elsevier: Amsterdam, 1987. (9) (a) Boerrigter, P. M.; te Velde, G.; Baerends, E. J. Int. J. Quantum Chem. 1988, 33, 87. (b) te Velde, G.; Baerends, E. J. J. Comput. Phys. 1992, 99, 84. (10) (a) Snijders, G. J.; Baerends, E. J.; Vemooijs, P. At. Nucl. Data. Tables 1982, 26, 483. (b) Vemooijs, P.; Snijders, G. J.; and Baerends, E. J. Slater Type Basis Functions f o r the whole Periodic System. Intemal report, Free University of Amsterdam, The Netherlands, 1981. (1 1) Krijn, J.; Baerends, E. J. Fitfunctions in the HFS-method. Intemal Report (in Dutch), Free University of Amsterdam, The Netherlands, 1984. (12) (a) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J . Phys. 1980, 58, 1200. (b) Becke, A. D. Phys. Rev. 1988, A38, 2398. (c) Perdew, J. P. Phys. Rev. 1986, B33, 8822; 1986, B34, 7046. (13) Fan, L.; Ziegler, T. J . Chem. Phys. 1991, 94, 6057. (14) Fan, L.; Ziegler, T. J . Chem. Phys. 1991, 95, 7401. (15) (a) CsBszL, P., Pulay, P. J. Mol. Struct. 1984, 114, 31. (b) Pulay, P. Chem. Phys. Lett. 1980, 73, 393. (c) Pulay, P. J . Compur. Chem. 1982, 3, 556. (16) (a) Fan, L.; Versluis, L.; Ziegler, T.; Baerends, E. J.; Ravenek, W. Int. J . Quantum Chem. 1988, S22, 173. (b) Fan, L.; Ziegler, T. J. Chem. Phys. 1992, 96, 9005. (c) Fan, L.; Ziegler, T. J . Phys. Chem. 1992, 96, 6937. (17) Maurer, F.; Wieser, H. The University of Calgary, Calgary, Canada. (18) Schachtschneider; Shell Development Co., Emeryville, CA, 1960. (19) Jost, A,; Rees, B.; Yelon, W. B. Acta Crystallogr. 1975, B31,2649. (20) Whitaker, A.; Jeffery, J. W. Acta Crystallogr. 1967, 23, 977. (21) Hedberg, L.; Lijima, T.; Hedberg, K. J . Chem. Phys. 1979, 70, 3224. (22) Beagley, B.; Schmidling, D.; J . Mol. Struct. 1974, 22, 466. (23) (a) Beagley, B.; Cruickshank, D. W. J.; Pinder, P. M.; Robiette, A. G.; Sheldrick, G. M. Acta Crystallogr. Sect. B 1969, 25, 737. (b) Alemnningen, A,; Haaland, A.; Whal, K. Acta Chem. S c a d . 1%9,23,2245. (24) (a) Stoicheff, B. P. Can. J. Phys. 1954,32,339. (b) Pliva, J.; Johns, J. W. C.; Goodman, L. J. Mol. Spectrosc. 1991, 148, 427. (25) Watson, J. K. G. J. Mol. Spectrosc. 1973, 48, 479. (26) Pulay, P.; Fogarasi, G.; Boggs, J. E. J . Chem. Phys. 1981,74,3999. (27) (a) Pongor, G.; Pulay, P.; Fogarasi, G.; Boggs, J. E. J. Am. Chem. Soc. 1984, 106, 2765. (b) Pongor, G.; Fogarasi, G.; Magdo, I.; Boggs, J. E.; Kereszturi, G.; Ignatyev, I. Spectrochim. Acta, Part A 1992, 48, 111. (28) Handy, N. C.; Maslen, P. E.; Amos, R. D.; Andrews, J. S.; Murray, C. W.; Laming, G. J. Chem. Phys. Lett. 1992, 197, 506. (29) Goodman, L.; Ozkabak, A. G.; Thakur, S. N. J . Phys. Chem. 1991, 95, 9044. (30) Jones, L. H.; McDowell, R. S.; Goldblatt, M. Inorg. Chem. 1969, 8, 2349. (31) Sosa, C.; Andzelm, J.; Elkin, B. E.; Wimmer, E.; Dobbs, K. D.; Dixon, D. A. J . Phys. Chem. 1992, 96, 6630. (32) Jones, L. H.; McDowell, R. S.; Goldblatt, M.; Swanson, B. I. J . Chem. Phys. 1972, 57, 2050. (33) (a) Brateman, P. S. Metal Carbonyl Spectra; Academic: London, 1975. (b) Van Rentergem, M.; Claeys, E. G.; Van Der Kelen, P. G. J . Mol, Struct. 1983, 99, 207. (34) (a) BBrces, A.; Maurer, F.; Wieser, H.; Ziegler, T. Unpublished work. (b) Fan, L.; Ziegler, J. Chem. Phys. 1991, 95, 7401. JP950853R