Mass scaling for vibrational frequencies from ab initio calculations

Chemical Physics Letters 403 (2005) 275–279 www.elsevier.com/locate/cplett

Mass scaling for vibrational frequencies from ab initio calculations Karl K. Irikura

*

Computational Chemistry Group 100, National Institute of Standards and Technology, Physical and Chemical Properties Division, United States Department of Commerce, Gaithersburg, MD 20899-8380, USA Received 23 August 2004; in final form 7 January 2005 Available online 20 January 2005

Abstract Ab initio calculations of harmonic vibrational frequencies deviate from observed vibrational spectra because of errors in the ab initio potential and anharmonicity. Based upon the Morse model potential, an effective mass can be used to mimic anharmonic effects. Exploratory calculations show that this strategy is effective for diatomic molecules. Parametric expressions are provided for predicting fundamental frequencies and zero-point energies from calculated harmonic frequencies. Ó 2005 Published by Elsevier B.V.

1. Introduction One of the most popular uses of quantum chemistry is to predict vibrational frequencies, which are used for spectroscopy, thermochemistry, and kinetics (zeropoint energies, partition functions, Wigner tunneling correction). A multiplicative Ôscaling factorÕ is frequently applied to the computed vibrational frequencies. This practice dates from the observation that harmonic frequencies from HF/3-21G calculations are about 12% higher than experimentally observed fundamental frequencies [1]. This corresponds to a scaling factor of s
0.89 in Eq. (1), where xcalc is the calculated harmonic frequency and mcalc is the predicted fundamental (i.e., anharmonic) frequency. Systematic studies have established empirical scaling factors for many computational models [2–5]. mcalc ¼ sxcalc :

ð1Þ

Mathematically, frequency scaling is equivalent to scaling the potential energy surface by s2 relative to the minimum. Better results can be obtained by using different scaling factors for different internal-coordinate *

Fax: +1 301 869 4020. E-mail address: [email protected].

0009-2614/$ - see front matter Ó 2005 Published by Elsevier B.V. doi:10.1016/j.cplett.2005.01.022

force constants, as championed by Pulay and co-workers [6–8]. Besides errors in the potential, the other principal error is the neglect of anharmonic effects. The scaling factor absorbs this to some extent. Adding a quadratic term to Eq. (1) improves the fit and has been rationalized in terms of anharmonicity [9]. In the present investigation, we consider using two adjustable parameters, one for scaling the potential (i.e., force constants) and one to correct for anharmonicity.

2. Theory Multiplicative scaling, as described above, is adopted here for correcting the shape of the potential. We call the parameter s in Eq. (1) the frequency-scaling parameter for consistency with the literature, but its physical significance is to scale the potential energy surface (by s2). To correct for anharmonicity, we consider the Morse potential as a model for bond stretching. It is given by Eq. (2), where r is the internuclear distance, re is the equilibrium bond length, De is the dissociation energy, and b is a parameter [10]. The corresponding energy levels are given by Eqs. (3)–(5) where xe and xexe are the harmonic frequency and the anharmonicity constant, respectively, and n is the quantum number [10]. In

276

K.K. Irikura / Chemical Physics Letters 403 (2005) 275–279

Eqs. (4) and (5), h is PlanckÕs constant, c is the speed of light, and l is the reduced mass. The fundamental vibrational frequency is m = G(1)  G(0) = xe  2xexe. Note that a small mass has a large amplitude of motion, accentuating anharmonic effects (Eq. (5)) 2

V ðrÞ ¼ De ½1  exp ðbðr  re ÞÞ ;
GðnÞ ¼ xe


2 1 1 nþ ;  xe xe n þ 2 2

y¼

 1=2 xe ¼ b2 De h=2p2 lc ;

ð4Þ

xe xe ¼ hb2 =8p2 lc:

ð5Þ

We seek an effective mass, m, which, when substituted for l in Eq. (4), will provide a fictitious harmonic frequency that coincides with the actual fundamental. That is, we choose m to satisfy Eq. (6). The appropriate value is given in Eq. (7), where the constant A has the value shown. Although A is not a simple multiplicative scaling factor, for simplicity of language we call it the mass-scaling parameter  2 1=2  2 1=2 b De h=2p2 mc ¼ b De h=2p2 lc  2hb2 =8p2 lc; ð6Þ

y min

.X  x2calc ; xcalc mobs X X x2calc : ¼ m2obs  s2

ð9Þ

ð10Þ

For mass scaling only, described by Eq. (8), the optimum scaling parameter and corresponding minimized sum are given by Eq. (11) P xcalc ðxcalc  mobs Þl1=2 P 2 1 ; A¼ xcalc l X X 2 y min ¼ x2calc l1 : ðxcalc  mobs Þ  A2 ð11Þ Both the potential and the mass can be scaled simultaneously, as given by Eq. (12). The corresponding parameter values are given in Eqs. (13) and (14), and the minimized sum is given in Eq. (15). For purposes of comparison, the root-mean-square ÔerrorÕ, RMSE = (ymin/N)1/2, is more convenient than the sum itself   mcalc ¼ sxcalc 1  Al1=2 ; ð12Þ P s¼

ð7Þ

The value of A depends upon the vibrational interval being mimicked. For example, to mimic the (n  1)th overtone, given by G(n)  G(0), A should be replaced by A(n + 1)/2 in the effective mass in Eq. (7). In particular, to mimic the zero-point energy, given by G(0), A should be replaced by A/4. The use of mass scaling is motivated by the observation for diatomic molecules that xexe µ xel1/2, approximately [11]. The form of Eq. (7) is derived from the Morse model potential. Actual stretching potentials deviate from this model, so A is considered an adjustable parameter here. The present, exploratory study is restricted to diatomic molecules. For a diatomic molecule, the correction can be applied to the harmonic frequency according to Eq. (8)   mcalc ¼ xcalc ðl=mÞ1=2 ¼ xcalc 1  Al1=2 : ð8Þ

ðmcalc;i  mobs;i Þ2 ;

X

s¼ ð3Þ

N X i¼1

ð2Þ



  m1=2 ¼ l= l1=2  A ;  1=2 ¼ 2l1=2 xe xe =xe : A ¼ b2 h=8p2 cDe

scaling factor and minimized sum are given in Eq. (10), where the summation index (i) and range (1 to N) are implicit

P  P 2 1  P x2calc l1=2 xcalc l xcalc mobs l1=2  ð xcalc mobs Þ ; P 2 2 P 2 P 2  xcalc xcalc l1 xcalc l1=2 

ð13Þ A¼ y min

P x2calc  xcalc mobs P 2 1=2 ; s xcalc l X ¼ mobs ðmobs  sxcalc Þ þ sA X  xcalc mobs l1=2 : s

P

ð14Þ

ð15Þ

It is not worthwhile to apply anharmonic corrections to crude calculations. Thus, the high-quality aug-cc-pVTZ basis sets [12,13] are used in the present study. The ccpVTZ basis is used for Li because no diffuse functions are available. CCSD(T) calculations are done using the ACES II software package [14,15]. All other computations are done using the GAUSSIAN package [15,16]. Core electrons (K-shell) are frozen in all post-Hartree– Fock calculations, except for Li centers. 4. Results and discussion

3. Computational details Parameter values are fitted by minimizing the sum of square differences between predicted (mcalc) and observed (mobs) fundamental frequencies for a set of N diatomic molecules, as described by Eq. (9). Simple frequency scaling is given by Eq. (1). The corresponding optimum

Variability of ÔconstantÕ. Mass scaling can only be useful if the parameter value is nearly constant across a representative group of molecules. To explore the variability of this parameter, it is computed for 16 diatomic molecules (HF, DF, HCl, DCl, H2, T2, Hþ 2 , HBr, CH, Cl2, CN, CO, CO+, CP, CS, Cs2) using Eq. (7) and

K.K. Irikura / Chemical Physics Letters 403 (2005) 275–279

experimental spectroscopic constants [17]. For this set, fundamental frequencies range from 42 to 4159 cm1. The resulting values of A range from 0.030 u1/2 (for CS) to 0.043 u1/2 (for HF and CH), with a mean value of 0.036 u1/2 and standard deviation of 0.004 u1/2. This is close to the mean values 0.040 u1/2 and 0.038 u1/2 previously obtained for some diatomic hydrides and nonhydrides, respectively [11]. In a set of 50 molecules from [11], only ClF and F2 show surprising values (A = 0.055 and 0.076 u1/2, respectively). Thus, we conclude that the parameter A is indeed reasonably uniform across diatomic molecules. Fundamental frequencies. The first test set consists of the six isotopic forms of dihydrogen: H2, HD, HT, D2, DT, and T2. As above, experimental spectroscopic constants are taken from the compilation by Huber and Herzberg [17]. Table 1 shows the frequency-scaling factors s (Eq. (1)) and RMSEs for full configuration interaction (FCI, performed as CCSD), the hybrid B3LYP density functional [18,19], second-order perturbation theory (MP2), and Hartree–Fock (HF). The Born– Oppenheimer potential curves are identical for isotopologues, so they require identical scaling factors in the harmonic limit. However, for frequency scaling, the RMSE is large, at 27 cm1. It is the same for all four theoretical potentials, confirming that it is due to anharmonicity. It has already been pointed out that frequency scaling is problematic for isotopologues [20]. The result of mass scaling alone (Eq. (8)) is also included in Table 1. The RMSEs for the FCI and B3LYP potentials are very good, indicating that anharmonic effects predominate. The values of the mass-scaling parameter A are also reasonable. When both scaling parameters are in-

277

cluded (Eq. (12)), all four theoretical potentials yield the same, excellent results (RMSE = 0.4 cm1) and the same reasonable value of the mass-scaling parameter (A = 0.040 u1/2). The values of A from the two-parameter fits are more reasonable than those obtained using mass scaling alone. The second test set consists of the 11 molecules FH, FD, FT, LiH, LiD, LiF, Li2, N2, CO, BeO, and O2. As shown in Table 2, frequency scaling works better than mass scaling except for the highest-quality potential, CCSD(T), and for MP2. In the two-parameter fits, the mass-scaling parameter has unexpected values except for the CCSD(T) potential. Fitting both parameters simultaneously is always at least as good as either parameter alone, as required mathematically, but is only dramatically better for the CCSD(T) and MP2 potentials. Results for this set of molecules, although good, are poorer than for the dihydrogen molecules, since they have diverse potential energy functions. Zero-point energies (ZPEs). The ZPE is typically taken as one-half of the sum of the computed (harmonic) vibrational frequencies. A multiplicative scaling factor is frequently applied to improve accuracy [2–4,21]. Table 3 shows results for frequency scaling, mass scaling, and their combination, for fitting the ZPEs of the dihydrogen molecules. The constants are given subscripts Ô0Õ to indicate that they are fitted to ZPEs instead of fundamentals. ZPEs are about half the magnitude of fundamentals, so the RMSEs are smaller. Since the anharmonic effects are eight times smaller than for fundamentals, the advantage of mass scaling over frequency scaling is diminished. In the twoparameter fits, the ratio A0/A = 0.28 is close to the

Table 1 Accuracy of predicted fundamental frequencies for dihydrogen isotopologues, based upon Eqs. (1), (8) and (12) (aug-cc-pVTZ basis set) Calc.

FCI B3LYP MP2 HF

No scaling

Frequency scaling

Mass scaling

RMSE

s

RMSE

A (u1/2)

RMSE

Two-parameter scaling s

A (u1/2)

RMSE

157.7 170.5 247.9 300.5

0.955 0.951 0.930 0.917

27.0 27.0 27.0 27.0

0.039 0.042 0.060 0.071

0.6 1.8 15.0 23.9

1.001 0.997 0.975 0.961

0.040 0.040 0.040 0.040

0.4 0.4 0.4 0.4

RMSE values in cm1.

Table 2 Accuracy of predicted fundamental frequencies for second set of diatomic molecules (aug-cc-pVTZ basis set) Calc.

CCSD(T) CCSD B3LYP MP2 HF

No scaling

Frequency scaling

Mass scaling

RMSE

s

RMSE

A (u1/2)

RMSE

s

A (u1/2)

RMSE

61.8 95.2 67.0 86.7 287.4

0.976 0.958 0.972 0.988 0.884

33.9 21.3 28.5 82.9 79.3

0.036 0.053 0.033 0.029 0.136

16.1 32.2 39.9 72.0 151.1

1.017 0.967 0.962 1.083 0.842

0.055 0.013 0.015 0.119 0.070

10.1 19.9 27.2 40.3 69.6

RMSE values in cm1.

Two-parameter scaling

278

K.K. Irikura / Chemical Physics Letters 403 (2005) 275–279

Table 3 Accuracy of predicted zero-point energies (ZPEs) for dihydrogen isotopologues (aug-cc-pVTZ basis set) Calc.

No scaling

FCI B3LYP MP2 HF

Frequency scaling

Mass scaling 1/2

Two-parameter scaling

RMSE

s0

RMSE

A0 (u )

RMSE

s0

A0 (u1/2)

RMSE

19.2 25.6 64.5 90.9

0.989 0.985 0.964 0.949

3.8 3.8 3.8 3.8

0.010 0.013 0.031 0.043

0.6 0.6 7.1 11.6

1.002 0.998 0.976 0.962

0.011 0.011 0.011 0.011

0.2 0.2 0.2 0.2

RMSE values in cm1.

Table 4 Accuracy of predicted zero-point energies (ZPEs) for second set of diatomic molecules (aug-cc-pVTZ basis set) Calc.

No scaling

CCSD(T) CCSD B3LYP MP2 HF

Frequency scaling

Mass scaling 1/2

Two-parameter scaling

RMSE

s0

RMSE

A0 (u )

RMSE

s0

A0 (u1/2)

RMSE

8.4 25.5 21.6 35.9 123.4

0.998 0.979 0.994 1.010 0.903

8.1 11.0 20.5 34.1 45.0

0.005 0.023 0.002 0.001 0.108

7.1 16.5 21.6 35.8 76.5

1.015 0.964 0.959 1.082 0.839

0.023 0.021 0.051 0.090 0.109

4.6 9.2 14.9 20.7 34.1

RMSE values in cm1.

value (1/4) expected for a Morse oscillator. However, the ratios range from 0.24 (for FCI) to 0.60 (for HF) when only the mass is scaled. Table 4 shows the results for ZPEs of the second set of molecules. As above, mass scaling is less advantageous for ZPEs than for fundamentals because of the smaller anharmonic effects. For mass scaling, the ratios A0/A now vary more widely, from 0.05 (for MP2) to 0.79 (for HF). For the two-parameter fit, the corresponding ratios vary greatly, from 1.6 (for CCSD) to 3.4 (for B3LYP). In all situations shown in Tables 1–4, it may be observed that A0
A  0.03 u1/2. (Curiously, this is the difference expected for a Morse oscillator with a typical value of A
0.04 u1/2.) This empirical relationship may be used, along with the approximation s0
s, to avoid fitting the ZPEs explicitly. This works well, as shown in Table 5. Comparison with Tables 3 and 4 shows that little is gained by fitting the ZPEs directly.

Table 5 Rms residuals (in cm1) for ZPE prediction using parameters estimated from fitting to fundamentals (aug-cc-pVTZ basis set), as s0 = s and A0 = A  0.03 u1/2 Theory

Mass scaling Set 2

Set 1

Set 2

1.1 1.3 7.6 12.1

7.1 16.5 21.6 35.9 76.5

1.0 1.0 1.0 1.0

4.6 9.3 15.0 20.7 34.3

a

CCSD(T) CCSD (FCIb) B3LYP MP2 HF a

Two-parameter

Set 1

There is no triples correction for set 1 because there are only two electrons. b CCSD is equivalent to FCI for set 1.

5. Conclusions For modest ab initio calculations of vibrational frequencies of diatomic molecules, the principal source of error is the potential energy function, so simple frequency scaling is recommended. For high-quality calculations, non-linear mass scaling is more effective for predicting vibrational fundamentals from computed harmonic frequencies, especially for isotopologues. Parameters appropriate for predicting anharmonic ZPEs can be estimated from the parameters for fundamental frequencies. Generalization to polyatomic molecules is under investigation. References [1] J.A. Pople, H.B. Schlegel, R. Krishnan, D.J. Defrees, J.S. Binkley, M.J. Frisch, R.A. Whiteside, R.F. Hout, W.J. Hehre, Int. J. Quantum Chem. Symp. 15 (1981) 269. [2] J.A. Pople, A.P. Scott, M.W. Wong, L. Radom, Israel J. Chem. 33 (1993) 345. [3] M.W. Wong, Chem. Phys. Lett. 256 (1996) 391. [4] A.P. Scott, L. Radom, J. Phys. Chem. 100 (1996) 16502. [5] M.D. Halls, J. Velkovski, H.B. Schlegel, Theor. Chem. Acc. 105 (2001) 413. [6] G. Rauhut, P. Pulay, J. Phys. Chem. 99 (1995) 3093. [7] G. Rauhut, P. Pulay, J. Phys. Chem. 99 (1995) 14572. [8] J. Baker, P. Pulay, J. Comput. Chem. 19 (1998) 1187. [9] H. Yoshida, K. Takeda, J. Okamura, A. Ehara, H. Matsuura, J. Phys. Chem. A 106 (2002) 3580. [10] G. Herzberg, second edn.Spectra of Diatomic Molecules, vol. 1, van Nostrand Reinhold, New York, 1950. [11] P. Hassanzadeh, K.K. Irikura, J. Comput. Chem. 19 (1998) 1315. [12] T.H. Dunning Jr., J. Chem. Phys. 90 (1989) 1007. [13] D.E. Woon, T.H. Dunning Jr., J. Chem. Phys. 98 (1993) 1358. [14] J.F. Stanton, et al., ACES II, release 3.0 Quantum Theory Project, University of Florida, 2002.

K.K. Irikura / Chemical Physics Letters 403 (2005) 275–279 [15] Certain commercial materials and equipment are identified in this Letter in order to specify procedures completely. In no case does such identification imply recommendation or endorsement by the National Institute of Standards and Technology, nor does it imply that the material or equipment identified is necessarily the best available for the purpose. [16] M.J. Frisch et al., GAUSSIAN 03 Version B.05, Gaussian Inc., Pittsburgh, PA, 2003. [17] Constants of Diatomic Molecules (data prepared by J.W. Gallagher and R.D. Johnson, III), K.P. Huber, G. Herzberg (Eds.),

[18] [19] [20] [21]

279

NIST Chemistry WebBook, National Institute of Standards and Technology, Gaithersburg, Maryland, 2001. Available from: . A.D. Becke, J. Chem. Phys. 98 (1993) 5648. P.J. Stephens, F.J. Devlin, C.F. Chabalowski, M.J. Frisch, J. Phys. Chem. 98 (1994) 11623. Y.N. Panchenko, V.I. Pupyshev, C.W. Bock, J. Mol. Struct. 550 (2000) 495. R.S. Grev, C.L. Janssen, H.F. Schaefer III, J. Chem. Phys. 95 (1991) 5128.