- Email: [email protected]
An approximation for hindered rotor state energies
Chemical Physics Letters 383 (2004) 203–207 www.elsevier.com/locate/cplett
An approximation for hindered rotor state energies John R. Barker
a,b,*
, Colleen N. Shovlin
a,b
a
b
Department of Atmospheric, Oceanic, and Space Sciences, University of Michigan, 1520 Space Research Building, 2455 Hayward Street, Ann Arbor, MI 48109-2143, USA Department of Chemistry, University of Michigan, 1520 Space Research Building, 2455 Hayward Street, Ann Arbor, MI 48109-2143, USA Received 24 October 2003; in final form 10 November 2003 Published online: 2 December 2003
Abstract A method suggested by Troe [J. Chem. Phys. 66 (1977) 4758] for estimating hindered internal rotation state energies is first evaluated and then extended by including anharmonic terms and introducing a smooth switching function. Comparisons between calculated thermodynamic functions and exact results [J. Chem. Phys. 10 (1942) 428] show that the extended method is greatly improved and performs at a level comparable to other approximate methods, but without their drawbacks. The extended method is very useful for calculating the densities of states and partition functions needed for reaction rate theories. Ó 2003 Elsevier B.V. All rights reserved.
1. Introduction Because hindered internal rotors (HIR) are ubiquitous, their theoretical description remains a subject of interest even 70 years after a quantum mechanical solution was first published [1]. Reviews and books (see citations by Truhlar [2], Knyazev [3,4], and McClurg et al. [5]) show that the problem has been solved essentially completely, and yet interest continues. This is because although the hindered rotor problem can be solved exactly, the solutions are cumbersome and timeconsuming. There is a need for methods that are more convenient and yet sufficiently accurate. Eigenstate energies are used for calculating sums and densities of states in statistical reaction rate theories [6–9] and for calculating partition functions. In 1991, Truhlar [2] described a simple approximation for calculating the partition function for a one-dimensional HIR, based on a multiplicative correction (the ÔPitzer– Gwinn ApproximationÕ [10]) to the harmonic oscillator partition function. From comparisons with the exact values from Pitzer and Gwinn (PG) [10], Truhlar found that the average relative error in the approximate par*
Corresponding author. Fax: +1-734-936-0503. E-mail address: [email protected] (J.R. Barker).
0009-2614/$ - see front matter Ó 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2003.11.024
tition function was 7%, which is acceptable for many applications. For applications that require sums and densities of states, not just partition functions, other methods must be used [3,5,4,11–14]. The purpose of this Letter is to report the accuracy of an extended version of a simple approximation recommended by Troe [15]. Troe noted the well-known fact [1,10,16] that the onedimensional HIR resembles a harmonic oscillator when the HIR energy is much lower than the potential energy barrier to internal rotation, and resembles a free rotor when the energy is much higher than the barrier. He recommended that the density of states be calculated as for a harmonic oscillator at energies below the barrier height and as for a free rotor at higher energies. Later Troe [17], used the classical phase space integral to develop semi-empirical approximations for classical sums and densities of HIR states. He also inverted the semi-empirical approximations to obtain approximate expressions for the quantum state energies. In a paper on unimolecular decomposition of the t-butyl free radical, Knyazev et al. [3] described an approximate method for calculating HIR densities of states, later corrected for quantum effects [4]. Essentially the same method was independently derived a few years
204
J.R. Barker, C.N. Shovlin / Chemical Physics Letters 383 (2004) 203–207
later by McClurg et al. [11], who were unaware of the previous derivation [12]. Partition functions calculated with this method are reasonably accurate: the maximum and average relative errors were reported to be <3% and <1%, respectively [4]. Although the partition functions are accurate, the estimated state energies are equally spaced, like a harmonic oscillator, but quite unlike hindered rotor states. In addition, an unphysical classical continuum appears. Subsequent publications discussed the choice of the zero of energy [13] and the degree of separability of HIR from over-all molecular rotation [4]. In the present work, we first examine the approximation recommended by Troe [15] in 1977 and show that it is only a little less accurate than TruhlarÕs method [2], and then we improve it. The present new method produces partition functions with accuracy comparable to the approximate methods of Knyazev [3,4] and McClurg et al. [5,11,14], but with much more realistic estimates of state energies. It is an advancement over TroeÕs approximations of quantum state energies [17] because its accuracy can be improved systematically by adding more terms. Sums and densities of states can be evaluated using the Beyer–Swinehart algorithm [18,19]. Ultimately, of course, exact numerical treatments, such as recent work by Van Speybroeck and colleagues [20– 23], will supplant all of the approximate methods, but until that time, the present method promises to be quite useful.
2. Theory 2.1. Hindered rotor The theory has been outlined in detail by many authors. Here, we follow the summary and much of the notation from McClurg et al. [5,14]. For a one-dimensional hindering potential given by a single Fourier component, the potential energy U can be written V ð1 cos n/Þ; ð1Þ 2 where / is the rotation angle, n is the number of minima, and V is the barrier height. When the rotor energy (E) is much smaller than V , the hindered rotor Schr€ odinger equation (equivalent to the Mathieu equation) reduces to that for a harmonic oscillator with frequency m: 1=2 n V m¼ ; ð2Þ 2p 2Ir
U¼
where Ir is the reduced moment of inertia for the HIR. The energy levels at E < V correspond to an anharmonic oscillator: ð3Þ Ev ¼ hm nv þ 12 þ Tv r; nv ;
where nv ¼ 0,1,2,3,. . . is the vibrational quantum number and Tv ðr; nv Þ is a sum of anharmonic terms [24,5]: ð1 þ 2nv þ 2n2v Þ þ Oðn3v =r2 Þ; ð4Þ 16r where r ¼ V =hm, and the sum of remaining terms is of the order of n3v =r2 . When the rotor energy E is much greater than the barrier height V , the hindered rotor Schr€ odinger equation reduces to that for a free internal rotor. The resulting state energies are those of a free rotor displaced in energy by V =2, plus additional terms: V Ej ¼ Bj2 þ þ Tr ðr; n; B; jÞ; ð5Þ 2 where j ¼ 0, 1, 2,. . ., is the rotational quantum number, B ¼ h2 =8p2 Ir is the rotational constant and Tr ðr; n; B; jÞ is the sum of additional terms of a series that is not uniformly convergent with increasing j. The terms in the series can be worked out from relationships involving the characteristic values of the Mathieu equation [5,24]. For example, when jjj > n=2, then Tv ðr; nv Þ ¼ hm
Tr ðr; n; B; jÞ ¼
r 4 n2 B 2
8½ð2j=nÞ 1
þ Oðr6 =j4 Þ:
ð6Þ
2.2. Approximation-1 Approximation-1 (suggested by Troe [15]) consists of assuming the state energies for E < V are given by Eq. (3) with Tv ðr; nv Þ ¼ 0, and state energies for E > V are given by Eq. (5) with Tr ðr; n; B; jÞ ¼ 0. Although not part of TroeÕs original description of this approximation, precautions must be taken to ensure that states are not inadvertently added or omitted in the vicinity of E V . This can occur because the state energies calculated from Eqs. (3) and (5) do not match. We have added such a precaution by numbering the quantum states (ns ¼ 1; 2; . . .) and obtaining quantum numbers nv and j in a manner that accounts correctly for degeneracies of the vibrational and rotational states: ns 1 nv ¼ INT ; ð7aÞ n n s j ¼ INT ; ð7bÞ 2 where INT(b) is the integer portion of positive real number b. We evaluated the performance of Approximation-1 by using quantum numbers obtained from Eqs. (7). The state energy Es was assumed to be given by Ev for Ev 6 V and Ej for Ev > V . For comparison with PG, the standard statistical mechanics formulae [10,25] were employed for free energy (F ), entropy (S), and molar heat capacity (C). Because volume is not involved in the thermodynamic functions for an internal rotation, Ôfree energyÕ (F ) can refer to either Helmholtz (A) or Gibbs
J.R. Barker, C.N. Shovlin / Chemical Physics Letters 383 (2004) 203–207
(G) free energy, and Ôheat capacityÕ (C) can refer to either constant volume (Cv ) or constant pressure (Cp ): F ¼ R ln Q; ð8aÞ T
Q1 S ¼ R ln Q þ ; ð8bÞ Q ( 2 ) Q2 Q1 ; ð8cÞ C¼R Q Q where R is the Gas Law constant and Q, Q1 , and Q2 are the hindered rotor partition function and its first two moments: 1 1X Es E0 Q¼ exp ; ð9aÞ n ns ¼1 kT 1 1X Es E 0 Es E0 Q1 ¼ exp ; n ns ¼1 kT kT
ð9bÞ
2 1 1X Es E 0 Es E0 exp Q2 ¼ : n ns ¼1 kT kT
ð9cÞ
In these equations, the zero of energy [13] is set equal to the energy of the lowest state (E0 , the zero point energy (ZPE)) for comparison with the PG tables. The factor of 1=n appears because of symmetry. The symmetry of each state depends in detail on the over-all angular momentum of the molecule and its projection on the molecular axis [26,27]. In general, only 1=n of the states have the proper symmetry to make the total wavefunction symmetrical. For a thermal distribution, a huge number of angular momentum states and projections are represented, making it impractical to keep track of them all. On average, all possible hindered rotor states are populated and thus we follow Waage and Rabinovitch [28] in calculating the partition function, its moments, and the sum of states by first utilizing all possible hindered rotor states and then dividing by n. This is the same strategy followed in dealing with free rotors, where the symmetry number is explicit. 2.3. Approximation-2 Approximation-1 is significantly improved by (a) including additional terms in the expressions for Ev and Ej and (b) incorporating a switching function to smoothly join the solutions for Es < V and Es > V . In the following, the first term in Tv ðr; nv Þ (Eq. (4)) is included in the solution for Ev and the first term in Tr ðr; n; B; jÞ (Eq. (6)) is included in the solution for Ej . The switching function SðEv ; V Þ is assumed to take the following form: SðEv ; V Þ ¼ ð1=2Þf1 þ tanh½aðEv V Þ=V g;
ð10aÞ
Es ¼ Ev ð1 SÞ þ Ej S:
ð10bÞ
205
By trial-and-error comparisons between the approximate and the exact state energies for methanol and ethane it was found that the results are not very sensitive to values of a in the range from a ¼ 4 to 6. We recommend a ¼ 5. This empirical function produces satisfactory results (see below), but other functions might prove to be even better choices. Further improvements can be achieved systematically by including more terms in the expressions for Tv ðr; nv Þ and Tr ðr; n; B; jÞ.
3. Results FO R T R A N subroutines incorporating Approximation-1 or Approximation-2 were used to calculate the state energies and the thermodynamic functions. The accuracy of thermodynamic functions calculated for a wide range of conditions is a general indicator of the performance of the approximation methods. In addition, comparisons are shown below between the approximate state energies and the exact values for methanol [27] and ethane [28]. PG constructed their tables by using two reduced variables [16,10]: V =RT and 1=Qf , where Qf is the partition function for a free rotor: 1 1=2 ½8p3 Ir kT : ð11Þ nh By using these variables, PR required only an accurate value for R, the Gas Law constant: all of the other fundamental constants are subsumed in the dimensionless variables. The Gas Law constant has undergone only slight revision since 1942: PG assumed R ¼ 1:9869 cal K1 mol1 , which agrees with the currently accepted value [29] (8.314472 J K1 mol1 ¼ 1.987207 cal K1 mol1 ) within 0.001 cal K1 mol1 . Therefore, no corrections were made in comparing the present results to the PG tables. The PG tables [10] cover the ranges 0 6 ðV =RT Þ 6 20 and 0 6 Q1 f 6 0:8, with some missing data in the region ðV =RT Þ 6 7 and Q1 f P 0:6. In the present work, comparisons were made to all of the data in the following PG tables: Table I (for F =T ), Table III (for S), and Table VI (for C). The thermodynamic results obtained using Approximation-1 and Approximation-2 are presented in Fig. 1 as differences from the PG tabulations (Error ¼ Xapprox XPG ). According to Eqs. (8) errors in F =T are due to errors in the approximate partition function (Q), while errors in S and C may also be due to errors in Q1 and/or Q2 . The figures show that, as expected, the approximations tend to approach the exact values in the limits: V =RT ! 0 and V =RT ! 1. The approximations are least accurate in the range 0:3 6 V =RT 6 3. Despite its simplicity, Approximation-1 is sufficiently accurate for many purposes. The maximum error in Qf ¼
206
J.R. Barker, C.N. Shovlin / Chemical Physics Letters 383 (2004) 203–207
Fig. 1. Errors (cal K1 mol1 ) in thermodynamic functions calculated using Approximation-1 and Approximation-2. The upper and lower limits of the data (vertical) axis in each plot correspond to the maximum positive and negative errors, respectively.
Both Approximation-1 and Approximation-2 consist of estimating the energies of the HIR states. Thus the sum of states is a Ôstair-stepÕ function with all steps of equal height, but with varying energy differences between states. This is illustrated in Figs. 2 and 3, where the sum of states is compared with values calculated from the exact hindered rotor state energies for methanol [27] and for ethane [28], respectively. (The sum of states indicates ÔfractionalÕ states because the sum is calculated using all possible hindered rotor states and
5 Koehler and Dennison Approximation-2
4 Sum of States
F =T is 0.13 cal K1 mol1 , corresponding to a maximum error of 6.8% in the partition function. The average error in the partition function is 2.1%. This is very similar to the performance of TruhlarÕs approximate method [2]. Errors in S and C are also relatively small: maximum errors of 0.17 and 0.27 cal K1 mol1 , respectively, and average errors of 0.07 and 0.1 cal K1 mol1 , respectively. This level of performance is acceptable for many applications because the relative contribution of a single hindered rotor to the total entropy or heat capacity is usually small and other assumptions and approximations may introduce errors larger than these. Approximation-2 performs significantly better than Approximation-1, with little additional computational effort. For all of the thermodynamic functions, the zone of V =RT values where significant errors occur is narrower for Approximation-2 than for Approximation-1. This improvement is particularly apparent at large values of V =RT . Typically, the maximum errors in the thermodynamic functions are reduced by factors of 2–3, while average errors are reduced by factors of 4–5, compared to Approximation-1. In particular, the maximum and average errors in F =T are 0.065 and 0.011 cal K1 mol1 , respectively, corresponding to errors in the partition function of 3.3% and 0.6%, respectively. This performance is similar to the Knyazev [3,4] and McClurg et al. [5,11,14] approximation, which has a stated maximum error of <3% [4].
3 2
zpe -1
V = 769.4 cm
1 0 0
500
1000 Energy (cm-1)
1500
Fig. 2. Comparison of the exact sum of states for methanol from Koehler and Dennison [27] to the results from Approximation-2. In both cases, all states are summed, regardless of symmetry, and the result is divided by the symmetry number (n ¼ 3). The ZPE and barrier height (V ) are shown.
J.R. Barker, C.N. Shovlin / Chemical Physics Letters 383 (2004) 203–207
207
Acknowledgements
Sum of States
8
Approximation-2 Exact (Waage & Rabinovitch)
We thank Lawrence L. Lohr for stimulating discussions. Thanks go to the Petroleum Research Fund (administered by the American Chemical Society) for partial financial support of this project. C.N.S. also thanks the College of Engineering at the University of Michigan for a Sarah Marian Parker Fellowship for undergraduate research.
6 4 zpe -1
2
V = 1024 cm
0 0
500
1000 Energy (cm-1)
1500
2000
Fig. 3. Comparison of the exact sum of states for ethane from Waage and Rabinovitch [28] to the results from Approximation-2. In both cases, all states are summed, regardless of symmetry, and the result is divided by the symmetry number (n ¼ 3). The ZPE and barrier height (V ) are shown.
then dividing by n, as described above.) It is clear that the agreement with the exact results is very good, but not of sufficient quality for spectroscopic purposes. The greatest error in state energy occurs when Es V , but Es asymptotically approaches the correct state energies when Es V . Because the errors are largest when Es V , the partition function and its moments are most strongly affected when V =RT 1, in agreement with the results shown in Fig. 1. Approximation-2 is suitable for calculating partition functions when 3% errors are acceptable. This level of error is of the same order, or smaller than errors introduced by common approximations adopted for the reduced moment of inertia (see East and Radom [30] for a quantitative assessment). Other common approximations include the assumption that the rotor is not coupled with other internal rotors [22] and vibrational modes, that it is not coupled with the over-all molecular rotation [4], and that the HIR potential energy can be represented with a single Fourier component. The method described here is easily calculated, predicts state energies with reasonable accuracy, and does not introduce non-physical effects like a classical continuum. It is easily employed in calculating sums and densities of states via the Beyer– Swinehart algorithm [18,19]. Subroutines based on Approximation-2 can be found in current versions of computer programs DenSum and Thermo: components of the free source MultiWell Program Suite (http://aoss.engin.umich.edu/multiwell/), which originally [31] did not include the HIR.
References [1] H.H. Nielsen, Phys. Rev. 40 (1932) 445. [2] D.G. Truhlar, J. Comput. Chem. 12 (1991) 266. [3] V.D. Knyazev, I.A. Dubinsky, I.R. Slagle, D. Gutman, J. Phys. Chem. 98 (1994) 5279. [4] V.D. Knyazev, J. Phys. Chem. A 102 (1998) 3916. [5] R.B. McClurg, R.C. Flagan, W.A. Goddard III, J. Chem. Phys. 106 (1997) 6675. [6] P.J. Robinson, K.A. Holbrook, Unimolecular Reactions, WileyInterscience, New York, 1972. [7] W. Forst, Theory of Unimolecular Reactions, Academic Press, New York, 1973. [8] R.G. Gilbert, S.C. Smith, Theory of Unimolecular and Recombination Reactions, Blackwell Scientific, Oxford, 1990. [9] K.A. Holbrook, M.J. Pilling, S.H. Robertson, Unimolecular Reactions, Wiley, Chichester, 1996. [10] K.S. Pitzer, W.D. Gwinn, J. Chem. Phys. 10 (1942) 428. [11] R.B. McClurg, J. Chem. Phys. 108 (1998) 1748. [12] V.D. Knyazev, J. Chem. Phys. 111 (1999) 7161. [13] R.B. McClurg, J. Chem. Phys. 111 (1999) 7163. [14] R.B. McClurg, J. Chem. Phys. 111 (1999) 7165. [15] J. Troe, J. Chem. Phys. 66 (1977) 4758. [16] K.S. Pitzer, J. Chem. Phys. 5 (1937) 469. [17] J. Troe, Chem. Phys. 190 (1995) 381. [18] T. Beyer, D.F. Swinehart, Comm. Assoc. Comput. Machines 16 (1973) 379. [19] S.E. Stein, B.S. Rabinovitch, J. Chem. Phys. 58 (1973) 2438. [20] V. Van Speybroeck, D. Van Neck, M. Waroquier, S. Wauters, M. Saeys, G.B. Marin, J. Phys. Chem. A 104 (2000) 10939. [21] V. Van Speybroeck, Y. Borremans, D. Van Neck, M. Waroquier, S. Wauters, M. Saeys, G.B. Marin, J. Phys. Chem. A 105 (2001) 7713. [22] V. Van Speybroeck, D. Van Neck, M. Waroquier, J. Phys. Chem. A 106 (2002) 8945. [23] P. Vanteenkiste, V. Van Speybroeck, G.B. Marin, M. Waroquier, J. Phys. Chem. A 107 (2003) 3139. [24] M. Abramowitz, I.A. Stegun (Eds.), Handbook of Mathematical Functions With Formulas, Graphs, and Mathematical Tables, U.S. Government printing Office, Washington DC, 1964. [25] D.A. McQuarrie, Statistical Mechanics, Harper & Row, New York, 1976. [26] B.L. Crawford Jr., J. Chem. Phys. 8 (1940) 273. [27] J.S. Koehler, D.M. Dennison, Phys. Rev. 57 (1940) 1006. [28] E.V. Waage, B.S. Rabinovitch, Int. J. Chem. Kinetics 3 (1971) 105. [29] P.J. Mohr, B.N. Taylor, J. Phys. Chem. Ref. Data 28 (1999) 1723. [30] A.L.L. East, L. Radom, J. Chem. Phys. 106 (1997) 6655. [31] J.R. Barker, Int. J. Chem. Kinetics 33 (2001) 232.
An approximation for hindered rotor state energies John R. Barker
a,b,*
, Colleen N. Shovlin
a,b
a
b
Department of Atmospheric, Oceanic, and Space Sciences, University of Michigan, 1520 Space Research Building, 2455 Hayward Street, Ann Arbor, MI 48109-2143, USA Department of Chemistry, University of Michigan, 1520 Space Research Building, 2455 Hayward Street, Ann Arbor, MI 48109-2143, USA Received 24 October 2003; in final form 10 November 2003 Published online: 2 December 2003
Abstract A method suggested by Troe [J. Chem. Phys. 66 (1977) 4758] for estimating hindered internal rotation state energies is first evaluated and then extended by including anharmonic terms and introducing a smooth switching function. Comparisons between calculated thermodynamic functions and exact results [J. Chem. Phys. 10 (1942) 428] show that the extended method is greatly improved and performs at a level comparable to other approximate methods, but without their drawbacks. The extended method is very useful for calculating the densities of states and partition functions needed for reaction rate theories. Ó 2003 Elsevier B.V. All rights reserved.
1. Introduction Because hindered internal rotors (HIR) are ubiquitous, their theoretical description remains a subject of interest even 70 years after a quantum mechanical solution was first published [1]. Reviews and books (see citations by Truhlar [2], Knyazev [3,4], and McClurg et al. [5]) show that the problem has been solved essentially completely, and yet interest continues. This is because although the hindered rotor problem can be solved exactly, the solutions are cumbersome and timeconsuming. There is a need for methods that are more convenient and yet sufficiently accurate. Eigenstate energies are used for calculating sums and densities of states in statistical reaction rate theories [6–9] and for calculating partition functions. In 1991, Truhlar [2] described a simple approximation for calculating the partition function for a one-dimensional HIR, based on a multiplicative correction (the ÔPitzer– Gwinn ApproximationÕ [10]) to the harmonic oscillator partition function. From comparisons with the exact values from Pitzer and Gwinn (PG) [10], Truhlar found that the average relative error in the approximate par*
Corresponding author. Fax: +1-734-936-0503. E-mail address: [email protected] (J.R. Barker).
0009-2614/$ - see front matter Ó 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2003.11.024
tition function was 7%, which is acceptable for many applications. For applications that require sums and densities of states, not just partition functions, other methods must be used [3,5,4,11–14]. The purpose of this Letter is to report the accuracy of an extended version of a simple approximation recommended by Troe [15]. Troe noted the well-known fact [1,10,16] that the onedimensional HIR resembles a harmonic oscillator when the HIR energy is much lower than the potential energy barrier to internal rotation, and resembles a free rotor when the energy is much higher than the barrier. He recommended that the density of states be calculated as for a harmonic oscillator at energies below the barrier height and as for a free rotor at higher energies. Later Troe [17], used the classical phase space integral to develop semi-empirical approximations for classical sums and densities of HIR states. He also inverted the semi-empirical approximations to obtain approximate expressions for the quantum state energies. In a paper on unimolecular decomposition of the t-butyl free radical, Knyazev et al. [3] described an approximate method for calculating HIR densities of states, later corrected for quantum effects [4]. Essentially the same method was independently derived a few years
204
J.R. Barker, C.N. Shovlin / Chemical Physics Letters 383 (2004) 203–207
later by McClurg et al. [11], who were unaware of the previous derivation [12]. Partition functions calculated with this method are reasonably accurate: the maximum and average relative errors were reported to be <3% and <1%, respectively [4]. Although the partition functions are accurate, the estimated state energies are equally spaced, like a harmonic oscillator, but quite unlike hindered rotor states. In addition, an unphysical classical continuum appears. Subsequent publications discussed the choice of the zero of energy [13] and the degree of separability of HIR from over-all molecular rotation [4]. In the present work, we first examine the approximation recommended by Troe [15] in 1977 and show that it is only a little less accurate than TruhlarÕs method [2], and then we improve it. The present new method produces partition functions with accuracy comparable to the approximate methods of Knyazev [3,4] and McClurg et al. [5,11,14], but with much more realistic estimates of state energies. It is an advancement over TroeÕs approximations of quantum state energies [17] because its accuracy can be improved systematically by adding more terms. Sums and densities of states can be evaluated using the Beyer–Swinehart algorithm [18,19]. Ultimately, of course, exact numerical treatments, such as recent work by Van Speybroeck and colleagues [20– 23], will supplant all of the approximate methods, but until that time, the present method promises to be quite useful.
2. Theory 2.1. Hindered rotor The theory has been outlined in detail by many authors. Here, we follow the summary and much of the notation from McClurg et al. [5,14]. For a one-dimensional hindering potential given by a single Fourier component, the potential energy U can be written V ð1 cos n/Þ; ð1Þ 2 where / is the rotation angle, n is the number of minima, and V is the barrier height. When the rotor energy (E) is much smaller than V , the hindered rotor Schr€ odinger equation (equivalent to the Mathieu equation) reduces to that for a harmonic oscillator with frequency m: 1=2 n V m¼ ; ð2Þ 2p 2Ir
U¼
where Ir is the reduced moment of inertia for the HIR. The energy levels at E < V correspond to an anharmonic oscillator: ð3Þ Ev ¼ hm nv þ 12 þ Tv r; nv ;
where nv ¼ 0,1,2,3,. . . is the vibrational quantum number and Tv ðr; nv Þ is a sum of anharmonic terms [24,5]: ð1 þ 2nv þ 2n2v Þ þ Oðn3v =r2 Þ; ð4Þ 16r where r ¼ V =hm, and the sum of remaining terms is of the order of n3v =r2 . When the rotor energy E is much greater than the barrier height V , the hindered rotor Schr€ odinger equation reduces to that for a free internal rotor. The resulting state energies are those of a free rotor displaced in energy by V =2, plus additional terms: V Ej ¼ Bj2 þ þ Tr ðr; n; B; jÞ; ð5Þ 2 where j ¼ 0, 1, 2,. . ., is the rotational quantum number, B ¼ h2 =8p2 Ir is the rotational constant and Tr ðr; n; B; jÞ is the sum of additional terms of a series that is not uniformly convergent with increasing j. The terms in the series can be worked out from relationships involving the characteristic values of the Mathieu equation [5,24]. For example, when jjj > n=2, then Tv ðr; nv Þ ¼ hm
Tr ðr; n; B; jÞ ¼
r 4 n2 B 2
8½ð2j=nÞ 1
þ Oðr6 =j4 Þ:
ð6Þ
2.2. Approximation-1 Approximation-1 (suggested by Troe [15]) consists of assuming the state energies for E < V are given by Eq. (3) with Tv ðr; nv Þ ¼ 0, and state energies for E > V are given by Eq. (5) with Tr ðr; n; B; jÞ ¼ 0. Although not part of TroeÕs original description of this approximation, precautions must be taken to ensure that states are not inadvertently added or omitted in the vicinity of E V . This can occur because the state energies calculated from Eqs. (3) and (5) do not match. We have added such a precaution by numbering the quantum states (ns ¼ 1; 2; . . .) and obtaining quantum numbers nv and j in a manner that accounts correctly for degeneracies of the vibrational and rotational states: ns 1 nv ¼ INT ; ð7aÞ n n s j ¼ INT ; ð7bÞ 2 where INT(b) is the integer portion of positive real number b. We evaluated the performance of Approximation-1 by using quantum numbers obtained from Eqs. (7). The state energy Es was assumed to be given by Ev for Ev 6 V and Ej for Ev > V . For comparison with PG, the standard statistical mechanics formulae [10,25] were employed for free energy (F ), entropy (S), and molar heat capacity (C). Because volume is not involved in the thermodynamic functions for an internal rotation, Ôfree energyÕ (F ) can refer to either Helmholtz (A) or Gibbs
J.R. Barker, C.N. Shovlin / Chemical Physics Letters 383 (2004) 203–207
(G) free energy, and Ôheat capacityÕ (C) can refer to either constant volume (Cv ) or constant pressure (Cp ): F ¼ R ln Q; ð8aÞ T
Q1 S ¼ R ln Q þ ; ð8bÞ Q ( 2 ) Q2 Q1 ; ð8cÞ C¼R Q Q where R is the Gas Law constant and Q, Q1 , and Q2 are the hindered rotor partition function and its first two moments: 1 1X Es E0 Q¼ exp ; ð9aÞ n ns ¼1 kT 1 1X Es E 0 Es E0 Q1 ¼ exp ; n ns ¼1 kT kT
ð9bÞ
2 1 1X Es E 0 Es E0 exp Q2 ¼ : n ns ¼1 kT kT
ð9cÞ
In these equations, the zero of energy [13] is set equal to the energy of the lowest state (E0 , the zero point energy (ZPE)) for comparison with the PG tables. The factor of 1=n appears because of symmetry. The symmetry of each state depends in detail on the over-all angular momentum of the molecule and its projection on the molecular axis [26,27]. In general, only 1=n of the states have the proper symmetry to make the total wavefunction symmetrical. For a thermal distribution, a huge number of angular momentum states and projections are represented, making it impractical to keep track of them all. On average, all possible hindered rotor states are populated and thus we follow Waage and Rabinovitch [28] in calculating the partition function, its moments, and the sum of states by first utilizing all possible hindered rotor states and then dividing by n. This is the same strategy followed in dealing with free rotors, where the symmetry number is explicit. 2.3. Approximation-2 Approximation-1 is significantly improved by (a) including additional terms in the expressions for Ev and Ej and (b) incorporating a switching function to smoothly join the solutions for Es < V and Es > V . In the following, the first term in Tv ðr; nv Þ (Eq. (4)) is included in the solution for Ev and the first term in Tr ðr; n; B; jÞ (Eq. (6)) is included in the solution for Ej . The switching function SðEv ; V Þ is assumed to take the following form: SðEv ; V Þ ¼ ð1=2Þf1 þ tanh½aðEv V Þ=V g;
ð10aÞ
Es ¼ Ev ð1 SÞ þ Ej S:
ð10bÞ
205
By trial-and-error comparisons between the approximate and the exact state energies for methanol and ethane it was found that the results are not very sensitive to values of a in the range from a ¼ 4 to 6. We recommend a ¼ 5. This empirical function produces satisfactory results (see below), but other functions might prove to be even better choices. Further improvements can be achieved systematically by including more terms in the expressions for Tv ðr; nv Þ and Tr ðr; n; B; jÞ.
3. Results FO R T R A N subroutines incorporating Approximation-1 or Approximation-2 were used to calculate the state energies and the thermodynamic functions. The accuracy of thermodynamic functions calculated for a wide range of conditions is a general indicator of the performance of the approximation methods. In addition, comparisons are shown below between the approximate state energies and the exact values for methanol [27] and ethane [28]. PG constructed their tables by using two reduced variables [16,10]: V =RT and 1=Qf , where Qf is the partition function for a free rotor: 1 1=2 ½8p3 Ir kT : ð11Þ nh By using these variables, PR required only an accurate value for R, the Gas Law constant: all of the other fundamental constants are subsumed in the dimensionless variables. The Gas Law constant has undergone only slight revision since 1942: PG assumed R ¼ 1:9869 cal K1 mol1 , which agrees with the currently accepted value [29] (8.314472 J K1 mol1 ¼ 1.987207 cal K1 mol1 ) within 0.001 cal K1 mol1 . Therefore, no corrections were made in comparing the present results to the PG tables. The PG tables [10] cover the ranges 0 6 ðV =RT Þ 6 20 and 0 6 Q1 f 6 0:8, with some missing data in the region ðV =RT Þ 6 7 and Q1 f P 0:6. In the present work, comparisons were made to all of the data in the following PG tables: Table I (for F =T ), Table III (for S), and Table VI (for C). The thermodynamic results obtained using Approximation-1 and Approximation-2 are presented in Fig. 1 as differences from the PG tabulations (Error ¼ Xapprox XPG ). According to Eqs. (8) errors in F =T are due to errors in the approximate partition function (Q), while errors in S and C may also be due to errors in Q1 and/or Q2 . The figures show that, as expected, the approximations tend to approach the exact values in the limits: V =RT ! 0 and V =RT ! 1. The approximations are least accurate in the range 0:3 6 V =RT 6 3. Despite its simplicity, Approximation-1 is sufficiently accurate for many purposes. The maximum error in Qf ¼
206
J.R. Barker, C.N. Shovlin / Chemical Physics Letters 383 (2004) 203–207
Fig. 1. Errors (cal K1 mol1 ) in thermodynamic functions calculated using Approximation-1 and Approximation-2. The upper and lower limits of the data (vertical) axis in each plot correspond to the maximum positive and negative errors, respectively.
Both Approximation-1 and Approximation-2 consist of estimating the energies of the HIR states. Thus the sum of states is a Ôstair-stepÕ function with all steps of equal height, but with varying energy differences between states. This is illustrated in Figs. 2 and 3, where the sum of states is compared with values calculated from the exact hindered rotor state energies for methanol [27] and for ethane [28], respectively. (The sum of states indicates ÔfractionalÕ states because the sum is calculated using all possible hindered rotor states and
5 Koehler and Dennison Approximation-2
4 Sum of States
F =T is 0.13 cal K1 mol1 , corresponding to a maximum error of 6.8% in the partition function. The average error in the partition function is 2.1%. This is very similar to the performance of TruhlarÕs approximate method [2]. Errors in S and C are also relatively small: maximum errors of 0.17 and 0.27 cal K1 mol1 , respectively, and average errors of 0.07 and 0.1 cal K1 mol1 , respectively. This level of performance is acceptable for many applications because the relative contribution of a single hindered rotor to the total entropy or heat capacity is usually small and other assumptions and approximations may introduce errors larger than these. Approximation-2 performs significantly better than Approximation-1, with little additional computational effort. For all of the thermodynamic functions, the zone of V =RT values where significant errors occur is narrower for Approximation-2 than for Approximation-1. This improvement is particularly apparent at large values of V =RT . Typically, the maximum errors in the thermodynamic functions are reduced by factors of 2–3, while average errors are reduced by factors of 4–5, compared to Approximation-1. In particular, the maximum and average errors in F =T are 0.065 and 0.011 cal K1 mol1 , respectively, corresponding to errors in the partition function of 3.3% and 0.6%, respectively. This performance is similar to the Knyazev [3,4] and McClurg et al. [5,11,14] approximation, which has a stated maximum error of <3% [4].
3 2
zpe -1
V = 769.4 cm
1 0 0
500
1000 Energy (cm-1)
1500
Fig. 2. Comparison of the exact sum of states for methanol from Koehler and Dennison [27] to the results from Approximation-2. In both cases, all states are summed, regardless of symmetry, and the result is divided by the symmetry number (n ¼ 3). The ZPE and barrier height (V ) are shown.
J.R. Barker, C.N. Shovlin / Chemical Physics Letters 383 (2004) 203–207
207
Acknowledgements
Sum of States
8
Approximation-2 Exact (Waage & Rabinovitch)
We thank Lawrence L. Lohr for stimulating discussions. Thanks go to the Petroleum Research Fund (administered by the American Chemical Society) for partial financial support of this project. C.N.S. also thanks the College of Engineering at the University of Michigan for a Sarah Marian Parker Fellowship for undergraduate research.
6 4 zpe -1
2
V = 1024 cm
0 0
500
1000 Energy (cm-1)
1500
2000
Fig. 3. Comparison of the exact sum of states for ethane from Waage and Rabinovitch [28] to the results from Approximation-2. In both cases, all states are summed, regardless of symmetry, and the result is divided by the symmetry number (n ¼ 3). The ZPE and barrier height (V ) are shown.
then dividing by n, as described above.) It is clear that the agreement with the exact results is very good, but not of sufficient quality for spectroscopic purposes. The greatest error in state energy occurs when Es V , but Es asymptotically approaches the correct state energies when Es V . Because the errors are largest when Es V , the partition function and its moments are most strongly affected when V =RT 1, in agreement with the results shown in Fig. 1. Approximation-2 is suitable for calculating partition functions when 3% errors are acceptable. This level of error is of the same order, or smaller than errors introduced by common approximations adopted for the reduced moment of inertia (see East and Radom [30] for a quantitative assessment). Other common approximations include the assumption that the rotor is not coupled with other internal rotors [22] and vibrational modes, that it is not coupled with the over-all molecular rotation [4], and that the HIR potential energy can be represented with a single Fourier component. The method described here is easily calculated, predicts state energies with reasonable accuracy, and does not introduce non-physical effects like a classical continuum. It is easily employed in calculating sums and densities of states via the Beyer– Swinehart algorithm [18,19]. Subroutines based on Approximation-2 can be found in current versions of computer programs DenSum and Thermo: components of the free source MultiWell Program Suite (http://aoss.engin.umich.edu/multiwell/), which originally [31] did not include the HIR.
References [1] H.H. Nielsen, Phys. Rev. 40 (1932) 445. [2] D.G. Truhlar, J. Comput. Chem. 12 (1991) 266. [3] V.D. Knyazev, I.A. Dubinsky, I.R. Slagle, D. Gutman, J. Phys. Chem. 98 (1994) 5279. [4] V.D. Knyazev, J. Phys. Chem. A 102 (1998) 3916. [5] R.B. McClurg, R.C. Flagan, W.A. Goddard III, J. Chem. Phys. 106 (1997) 6675. [6] P.J. Robinson, K.A. Holbrook, Unimolecular Reactions, WileyInterscience, New York, 1972. [7] W. Forst, Theory of Unimolecular Reactions, Academic Press, New York, 1973. [8] R.G. Gilbert, S.C. Smith, Theory of Unimolecular and Recombination Reactions, Blackwell Scientific, Oxford, 1990. [9] K.A. Holbrook, M.J. Pilling, S.H. Robertson, Unimolecular Reactions, Wiley, Chichester, 1996. [10] K.S. Pitzer, W.D. Gwinn, J. Chem. Phys. 10 (1942) 428. [11] R.B. McClurg, J. Chem. Phys. 108 (1998) 1748. [12] V.D. Knyazev, J. Chem. Phys. 111 (1999) 7161. [13] R.B. McClurg, J. Chem. Phys. 111 (1999) 7163. [14] R.B. McClurg, J. Chem. Phys. 111 (1999) 7165. [15] J. Troe, J. Chem. Phys. 66 (1977) 4758. [16] K.S. Pitzer, J. Chem. Phys. 5 (1937) 469. [17] J. Troe, Chem. Phys. 190 (1995) 381. [18] T. Beyer, D.F. Swinehart, Comm. Assoc. Comput. Machines 16 (1973) 379. [19] S.E. Stein, B.S. Rabinovitch, J. Chem. Phys. 58 (1973) 2438. [20] V. Van Speybroeck, D. Van Neck, M. Waroquier, S. Wauters, M. Saeys, G.B. Marin, J. Phys. Chem. A 104 (2000) 10939. [21] V. Van Speybroeck, Y. Borremans, D. Van Neck, M. Waroquier, S. Wauters, M. Saeys, G.B. Marin, J. Phys. Chem. A 105 (2001) 7713. [22] V. Van Speybroeck, D. Van Neck, M. Waroquier, J. Phys. Chem. A 106 (2002) 8945. [23] P. Vanteenkiste, V. Van Speybroeck, G.B. Marin, M. Waroquier, J. Phys. Chem. A 107 (2003) 3139. [24] M. Abramowitz, I.A. Stegun (Eds.), Handbook of Mathematical Functions With Formulas, Graphs, and Mathematical Tables, U.S. Government printing Office, Washington DC, 1964. [25] D.A. McQuarrie, Statistical Mechanics, Harper & Row, New York, 1976. [26] B.L. Crawford Jr., J. Chem. Phys. 8 (1940) 273. [27] J.S. Koehler, D.M. Dennison, Phys. Rev. 57 (1940) 1006. [28] E.V. Waage, B.S. Rabinovitch, Int. J. Chem. Kinetics 3 (1971) 105. [29] P.J. Mohr, B.N. Taylor, J. Phys. Chem. Ref. Data 28 (1999) 1723. [30] A.L.L. East, L. Radom, J. Chem. Phys. 106 (1997) 6655. [31] J.R. Barker, Int. J. Chem. Kinetics 33 (2001) 232.