- Email: [email protected]
Upper and lower bounds to partition functions for molecules with double minimum potentials
.,. ._
.
:
. . .-‘.‘ UPPER AND LOWER ,ODNDS TO PARTITION FUNCTIONS -..: FGR MOLECULES : WITH-DOUBLE MINIMUM POTENTIALS .i._ __ .- .I. GRONHOLZ and W. WITSCHEL Departmht df C’henii.r&, UniirersiIy of Ulm. Germany .: Rx&d
i2 September 1977
Up&r and lower bounds to the viirational partition function 4 for molecules with double minimum potentials are derived. Both the exact lower (Gibbs-Bogoliubov) and upper (Golden-Thompsck) bound to q are evaluated analytically for the harmonic oscillator perturbed symmetrically or asymmetrically by a gaussian barrier. For the quadratic-quartic oscillator only the lower bound can be evaluated analytically, whereas the upper bound leads to a strongly divergent series resisting conventional summability techniques. In this case the classical partition function is used as an upper bound.
I. Introduction In a recent review on statistical methods for calculating thermodynamic functions, examples were given for molecules with double minimum vibrational potentials .[l] . Similar potentials are also interesting for hydrogen bonding and proton transfer [2,3] and for structural phase transitions in solids [4]. The customary technique is the separation of the remaining vibrational degrees of freedom and the numerical solution of the one- or two-dimensional Schriidinger equation either by direct solution
of the differential
tion [S] or by matrix dlagonallzation.
equa-
The molecular partition function q can then be evaluated by summation of the Boltzmann-factors calculated from the numerical eigenvalues. From a theoretical point of view it is more satisfactory to get analytical results which can be discussed generally. For a conventional force field the partition function can be obtained by means of thermodynamic perturbation theory leading in each order to simple results in ternis of exponential and hyperbolic functions [6]. As perturbations which cause the double minimum character of the potentialcannot be considered as small, different non-perturbatlonal.meth-odsmust beused. -Variational methods are non-p&turb&ive and have found wide.application in statistical mechanics. The recent reviews, by.Girardeau and Mazo [7] and
by Falk [8] summarize the mathematical foundation and applications of upper and lower bounds. Feynman [9] applied the lower bound to a linear oscillator with cubic anharmonicity which seems to be the only application to molecular vibrational problems. In this article we derive upper and lower bounds to the partition function q for the double minimum potentials mentioned. We shall not discuss again the mathematical background especially of the GoldenThompson inequality but refer to recent publications on this topic. The potentials behave well so that no principal difficulties arise. Therefore, only a brief sketch of the inequalities and the formalism to be used will be given in section 2. Section 3 contains the derivation of the bounds. The results are.evaluated numerically and discussed in section 4. The paper closes with a conclusion in section 5.
2. Variational principles and computational
methods
As the derivation of the following principles can be found in the literature [7-l I], we omit a repetition and just state the Gibbs-Bogoliubov and the Golden-Thompson inequality for the lower and upper bound to the partition function q. Consider a hamiltonian H, iiihich can be separated into
-’: .. _-: :
where fit, is not necessarily small. The Gibbs-Bogoliubov inequality represents.an upper.bound for the Helmholtz free energy A A=-~-llnq&40+(H1+r,
(2)
where 0-l = kT. A0 = -/1-l In-q0 is the_HehnhoItz free energy corresponding to HO, and (HI jn means - Tr{fil exp (-pfio)] _ @I& = (3) Tdexp t--P&>> Taking exponentials on both sides of eq. (2) yields the fohowing inequality for the partition functions 4 and qo q ~Tr{exp(-/I@)}
>q,-, exp (-@$~,)
=qL,_
(4)
which is the starting point for the calculations of lower bounds to ~7(the index “L” signifies “lower bound”). The Golden-Thompson inequality is 4 G qo(exp (-PHI PO = qu,
&xp[i(&
+ yj)] >n = exp (-$o~~ci>exp i-+y2y),
(10)
from which we .ob tain (11)
with (5)
(6)
In addition to the reviews mentioned above, we refer to the work of Breitenecker and Griimm [12] and of Ruskai [13] , dealing with the mathematical foundation of the Golden-Thompson inequality. Combining the inequalities, one has 4L =Q=Qu.
so that *r contains tberemainmg potential-energy terms.- This choicehas. the $dkta&; :&tit the &IcuIation of the &eragesm.@r>Oand
(exp (cZ))~ = exp <$fx2r-(>.
which will be used to compute upper bounds to 4. Taking Iogarithms on both sides one has a lower bound to the Heimholtz free energy A A >A, - fl-‘(exp (-&YL)J,.
’
(7)
The separation oftiinto two parts is fairly arbitrary, but from a practical point of view,.the number of possibilities is limited by the methods for evaluating. the traces in a suitable basis. If, for_example, one takes the kinetic energy operator for Ho and the potential energy operator for HII the averages can be calculated in the basis of plane waves with appropriate boundary conditions. Golden [lo] identified this qu as the cIassicaI partition function, qd- Therefore, one has another useful inequality
which wih be used to calculate 4” for the quartic OScihator. In ah other cases treated here, the appropriate
p = coth($&).
(12)
In the following we abbreviate /9zO by e. Formula (11) is the master formula for the subsequent calculations. It will be applied to the averages for the special cases we are interested in.
3. Calculation of the thermal averages 3.1. Potentials with gaussian barriers T&o types will be treated: ti‘asymmetric minimum potential with hamiltonian (I)
double
fi=fi0+&X2exp[-y(Z-a)2],
(I31
finding application in the theory of the hydrogen bond [2,3], and a symmetric double minbum potential with hamiltonian. (11) fl=Ho
+.(Yh~exp(-#),
.’ . . .
(14).
which is standard in the. molecular speckoscopi- of 1 large amplitude motions 116,171 .-Prelimbiaj c&u-lations using (II).as a model system in structtirz&ph&e transitions are prom&g. The principal forni of &se potentials is given in figs::l:and 2 bikmaji-v&y cotisiderably depending on the ~dten&d parameters.’ .. The thermal averages corres@%&ng to.(I) are given
287
.I. Gronholz, W. Witschel/Upper and lotier bounds to the molecularpartition function
power ofi appears in the exponential, and our master formula (11) cannot be applied directly. This difliculty is solved using the F,!urier-transformation exp[-r&i - A)2] = (4nn y)-lj2 +Ol X exp (--IrA-?/4ny) exp (i&) d7, s --m
(17)
which is a special applickion of Weyl’s prescription for quantizing classical functions, discussed extensively by Wilcox [18]. Formula (17) is introduced in eq. (16). The order of integration and trace formation is interchanged and Bloch’s theorem is applied -2 -1 COORDINATE
i
(clim&ionl&l
Fig. 1. Potential enera V(X) of hamiltonian (I). V(X) is the .sum of the harmonic oscillator potential V,(X) = (r!n/2)@ and the gaussian barrier V,(X) = fi 0 CLexp [ -T(X-A)~ ]_ Potential parametersare: (L= 2.7 = 1, A = 0.5.
(exp [-KY@? - A)‘] lo = (4 71~y)-l12 +m X exp (-h-A--?/4rry)(exp (iTJ@, dr s -m = (4 7111 y)-lj2
by
= (1 +izyp)- “’ :;
Gljo = k~iX’2exp [-r(_%! -
A)2] j,,
(15)
and (exp (-@,)+,
7
exp (-irA-T2/4rzy
.
- T2p/4) dr
[-A” /zy/( 1+izyr_1)]_
(18)
Combining this formula with eqs. (4) and (5); the upper and lower bounds to 4 for case I are
.
= (exp {-m exp [-~(2
4:’= q. expC--ad1+w)-~‘~exp[-A2d( 1+w)l 1, (1%
- A)2] }jo
=n$oy(exp[-rry(2-A)2])o,
(16)
*$I
=qo
5 I29
n-0
(1
t.yfi)-1’2
n!
X exp [-A2n-yl( 1 + JZYP)].
so that we are left with the calculation of (exp[-ny(* - A)2] jO_Unfortunately, the second
(20)
The corresponding formulae for the symmetric case II are obtained as special cases by putting A = 0. Thus, we get 4y
= q.
($0 = q.
ewI-4+rd-1121,
(21)
5 k$n.
(72)
( ltay/L)-“2.
n-0
The series for the upper bounds are convergent for any e, OLand 7 (7 > 0). The situation is less satisfactory for the quartic osci!lator as will be seen in section 3.2. COORDl&fATE
X (dimensionless1
Fig. 2. same as fig. 1, except that LI = 0 [hamiltonian (II)] in thgcase.
3.2. Quadratic-quartic
oscillator
The quadratic-quartic oscillator is of wide interest in quantum field theory, solid state physics, and statis-
288
J. Gronf~olz, W. .Wiisc?ze!lUpperand rower bounds to the moIecu~a~~artiti&fanHim
.tical mechanics. In molecular physics it is used to describe the strongly anharmonic puckering vibration of small ring molecules like cyclobutane [ 19] and oxetane [20,21] _Its hamiltonian is H=ti0m2(f&f4~).
(23
If& = -f , we have a pure quartic potential (case III) and iff? = 0, we simply add an z4-term to the harmanic oscillator potential (case IV). The averagesu and (exp(-/3Gr))u are ui, )u = tr L! *
(24)
ff4G?>,),
and (exp (&I)),
= (exp[--ECf2J?2 tf,J?)]
lo
:
-~
.-.
Instead of expanding exp (-&I), one. may try to transform exp (-@gl) into anotherexpression, which can tie.averaged by means-of ouimasier formula; Using the Weyl prescription,is.given by Wilcox .[18], ._ -. we have to cair$late the Fourier integral g(~)
=$. J
exp I-eV,X2
’
+&x4) exp (-irx)dx.
-a.
: (30) However, the Fourier transform of the function exp (-X4) could not be found in the tables of Fou:. rier transforms [23] and, theiefore, this m&hod fails, too. From the above results it is clear, that we cannot average the operator function exp[-eCf2p2 +f,J!‘!)] using Bloch’s theorem. As the classical partition function qcl was shown to be an upper bound to q, we
calculate 1 %l=G (W
The expressions dm>u are deduced from our master formula, if we‘expand the exponentials and compare equal powers of (Yon both sides of (11). We get
(26)
i-a
ss
exp (-fed)
--m
X expf_-_4e[(1+2f2)X2 +2f,X4]}dPdX,
~cl=(~7r/e)1’2[(i X e’ K1,4(4
Inserting this into eqs. (4) and (5), we have qL = q. exp ]--E($ f2c1 f %f,&l
,
(27)
and
(31)
which can be separated into a translational and a configurational part. Both integrals are available from an integral table [24] . The result is +2f2)/2f4]112 -
Q4W1,
(32)
where z is equalto ~(1 + 2f2)*/(32f4). The 1, are modified Bessel functions of first kind and fractional order V. In the case of a pure quartic oscillator, the above expression simplifies to &II) = aI12 2-r f,-‘/4&~(3/4)-~_
x
(4n - 2/c)!
2n-k
(2n - k)! 42n-k ’
-
(28)
(33)
Eq. (32) is not new. Using confluent hypergeometric functions instead of the I,, it was applied by Wehner and Bear-My1 [25] to calculate the zero field susceptibility of a system of classical anharmonic oscillators.
For case IV, qr, reduces to
(4’*Y
3.3. Improvement.of
(2%
>z!(~,z)!42n Unfortunately, the inftite series in eqs. (28) and (29) diverge-for any E, f2 and f4 (E,f4 # 0). Customary summation methods, which were successful for the Rayleigh-SchrEdinger perturbation series of the quartic oscillator [22] , fail as the series are too strongly divergent.
the botinds
I Our bounds can be improved by the.introduction of variational parameters using,a suggestion of Feyn-
man [9]. . The hamiltonian is separated into fir, and tiI_ &=r-I, +fil_
(34)
A potential energy term TA =h a’(.? - A’)2 is added
J. Gronhols~W. Witschel/Upper and lower bounds to the molecular partition function
and subtracted (35) and-the aver&cs~(fjl - F$) and (exp [-/3(kl - &)]) are performed, using the density.operator corresponding to (fro f fh)_ The resulting expressions can be
Table 1 C&es I and IL Dependence of the percent difference between upper and lower bound, 070, on the potential parameters (Y, ‘y,and A. D% is related to the arithmetic mean $7” + qL) and given by eq. (37). fin/W = 1
varied under the conditions 3qll(cr)/&
a0
(a = !2’or A’; 7 = U or L.).
(36)
For symmetric potentials one finds, as expected, Abptimum = 0. In all other cases considered here, the extremum conditions (36) lead to transcendental equations which cannot be solved analytically, and
therefore, the variation has to be performed numerically. This is ti disadvantage, as we are mainly interested in closed results. Numerical work is postponed to ;1 later paper, dealing with coupled osciilators. Only the qtiartic oscillator will be given as an example, because in this case the improvements from the variational procedure are most significant_ The expression to be optimized with respect to CL’is then given by (27) if we putf2 = -l/2 and C2= .Q’.
4. Numerical
results
4. I. Po terrtials con taikng gatcssian bakers
In tables 1 and 2 the difference between upper and lower bound, D = q u - qL, is used as a measure for the accuracy with which we can estimate 4. combining qu and qL. Related to the arithmetic mean f(qu+qL), we express D in percent D% = 2OOD/(q,
+ qL).
289
2 2
2 6
0 0
2 6
10 2 6 10 2 6 10 2 2 2 10
0 0 0 0 0 0 0 2 4 2 4
-6 6 10 10 10 2 2 10 10
11.88 9.31 7.79 63.42 47.58 39.59 90.86
0.4034 0.5621 0.6302 0.0713 0.1929 0.2719 0.0126 0.0667 0.1173 0.6572 0.9298 0.1447 0.3506
76.58 66.91 7.12 0.62 62.31 38.12
0.5122 0.6774 0.7368
0.3186 0.5430 0.6283 0.2632 0.4988 0.5917 0.7580 0.9413 0.6231 0.7827
bounds are best for small e-values, i.e. high tempera-
tures and/or low vibrational frequencies. Finally, in table 3,qU and qL are compared to a partition function 4 calculated by summation of Boltzmann-facts,rs with numerical eigenvaules from the literature. The comparison shows that qu is always cIoser to CJthan qL. This means that the upper bound is a better approximation to 9 than the lower bound. 4.2. Quark-quadratic
potettriol
Though our results apply for any combination
of
f2and f4 provided f4> 0, we confine our numerical
(37) Table 2.
Table 1 presents LEG and the upper and lower bound to the partition function 4 for cases I and II. E is kept fEed (E = 1) and &potential parameters (Y,7 and A
CaseI.
are varied. One sees, that II% decreases with decreasing CYand increasing y and A. Small a; large y and A stand for a “small” I?l in the sense of ordinary perturbation theory. This means, our bounds improve with decreasing magnitude of the perturbation. The next question is, wh>ther qu and qL approx@ate q better for high or for low temperatures. The
paramctcrs:
answer can be seen from table 2, where the relative difference of the bounds is listed again for various evalues and a given perturbation. It appears that our
Dependence of the percent difference D% between upper and iower bound on E = ML//CT. D% is related to the arithmetic mean &J + 4~) and given by eq. (37). Potential OL= Y = A = 2
(0
E
QL
D%
Q{f)
0.2
4.646 2.107 1.281
0.4 1.5 3.2 5.1 7.1 16.2 33.0
4.680 2.170 1.364 0.980 0.758 0.329 0.104
0.4 0.6
0.8
0.884
1.0 2.0 4.0
0.657 0.237 0.053
:.
290
J: Gronholz, W. iVitschei&pper and Iower bounds to t&molecular
: .’ : -.. .’ f :. I_:. .. -:_.. .:
p&titio~ function-,
..
Table 3
-.
L- ,.
--‘_
-_:
COIIIPirkm of the upper and lower bound to the partition function q for cases I and IL q is calculated by summmg~u~ the
(I, 11) qL
A
Y
0.289
19.05
0
0.578 0.868 5.844
9.53 6.35
0 0
29.440 29.440
0.025 0.025
29.440
0.025
0.159
0.9177 0.8473 0.7655 0.6211 0.4123 0.3916 0.3809
0 -0.505 -0.429 -0.382
x x x x
:.
1c+ lo-l2 1O-‘2 lo-‘2
.. ,(r,IIj
~.
0.9190 0.8548 0.7845 0.1137 0.8136
a) a) a) I x 10-t x lo-‘t .0.6235.x lo-” 0.5333 x IO-”
:.- (I,II) 4rJ ._ :
... :
.‘-:. 0.9209. 0.8614
:. -:
.I.
. b) b) b) b1
0.7969. : 0.1266 x 10-t 0.1440 x lo-to
.-0.1095 x io-‘0 ‘0.9469 x lo-”
a) Eigenvalueswere taken from Coon et al. [16]. b) Eigenvalucs were taken from Flanigan and de la Vega [2] -
investigation to the quartic oscillator because its eigenvalues are known best [26,27] and the computation of the modified Bessel functions [l/4 in (32) is superfluous. In table 4 we list qrD and the corresponding “variationaIly improved bound”, called qzft3. The upper bound is given by the classical partition function qcl. qcl approaches 4 with decreasing E. This is exactly the high-temperature behavlour to be expected for a classical partition function_ In contrast to qcl, the lower bounds approximate 4 best in a medium e-range for E = 1. The variations
from case I, II and also IV (though not treated here) and arises from the fact that the hamiltonian (III) does not contain any harmonic contribution to the potential energy. The best result for &tl deviates approximately 3% from the numerical 4.. 4.3. Summary We summarize the results. In cases I and II our bounds are quite accurate, if the height of the gaussian barrier is smaller than or comparable to the.Ievelspacing of the corresponding harmonic oscillator (i.e. (Y< 1) and if E < 1. For larger perturbations and evalues the results are still useful but less accurate_ For a purely anharmonic potential as m case III, we cannot abandon the numerical variation, because the lower bound without variation gives very poor. results. On the other hand, it is advantageous, that qLv gives good estimates of 4 even. for extremely strong anharmonicity by,a minimum of computational effort. The numerical variation is very simple and can be performed-even on8 programmable pocket calculator. The upper bound, q&s only useful for. high temperatures, i.e. for E Q 1. -_
lead to considerable improvements, especially for high temperatures, where qLv and qL differ by a factor
of IO*-103. The lower bound without variation is nearly useless in this case. This situation is different Table 4 Comparison of upper and lower bounds to the partition function of the quartic oscillator, q(I1ll. &ll) and #1 arc givenby (27) and (33). E=trR/kT. Forsfff), see text l
(III) qL
0.2
0.004
0.4 0.6 0.8 1.0
0.877 0x939 0.2579
2.0 4.0 6.0 8.0
0.2847 0.2251 0.0774 0.0232 0.0067
(III)
q(III) al
,(fU c-_
2.7182 1.5892 1.1469 0.9002 0.7386 0.3580 0.1153 0.0390 0.0132
2.8450 1.6593 1.1943 0.9350 0.7654 0.3693 0.1204 0.0416 0.0144
2.8757 1.7099 1.2616 1.0167 0.8600 0.5113 0.3041 0.2243 0.1808
4LV
5. Conclusion--
:
:
1.
:: We applied variational techniques to. the partition: functiodqf one-diniensionid strongly-anh&mon& model systems, finding application-iii molecular physics, but also m-other.fieIds__: ‘: 1. :-.~’.I.,. -: .I- .- ‘.. : . The present work’airnkd $t two pomts: .First,: to :. .~
a) Calculated from eigenvahtes of Bell et al. [27] _
I
.:__
. .
i
._
-,_
_:
--..
_
_I.G~o&olz. W. WitscheljUpper and iower bounds to the krolecutar partition finction
formulae ‘for the one-dimen&I& double minirnUni problems, second, to test the appl&bi& of variational methods t6. the partition function of coupled tihamlonic 6scillators. The calculation of eigetivalues and eigenfunctions [28] and the general b&aviour [29] of such systems is a field of active reSearch with many open mathematical questions_ For the solution of chemical problems, especially in the field of reaction kineitcs, it is necessary to find good approximations to the partition function. _The satisfying results for one-dimensional strongly anharmonic systems encourage calculations for coupled oscillators.
‘derive u&Xul anal&l
Acknowledgement Constructive remarks on an earlier version of this paper by Professor R.A. Sack are gratefully acknowledged. We also wish to thank the “Fonds der Chemischen Industrie” for financial support.
References [I ] S.G. Frankiss and J.H.S. Green, in: A Specialist Period[21 131 i4j [S] [6]
ical Report, Chemical Thermodynamics, Vol. 1 (The Chemical Society, London, 1973). M.C. Flanigan and J.R. de la Vega, J. Chem. Phys. 61 (1974) 1882. R.W. Redding, J. Mol. Spectry. 62 (1976) 8. H. Beck, J. Phys. C9 (1467) 33. J.W. Cooley, Math. Comput. 15 (1961) 363. W. Witschel, 6. Grosswendt and J. Bohmann, PTBMitteihmgen 87 (1972) 3.
291
[7] M.D. Girardcau and R.M. Mazo, Advan. Chcm. Phys. 24 (1973) 187. [8] H. Faik, Am. J. Phys. 38 (1970) 858. [9] R.P. Feynman, Statistical mechanics (Benjamin, New York, 1972). [lo] S. Golden, Phys. Rev. B 137 (1965) 1127. [ll] C.J. Thompson, J. Math. Phys. 6 (1965) 1812. [12] M. Breitcnecker and H.R. Griimm, Commun. Math. Phys. 26 (1972) 276. 1131 M.B. Ruskai, Commun. Math. Phys. 26 (1972) 280.
[14] F. Bloch, Z. Phgsik 74 (1932) 295. [lS] A. Messiah, Quantum mechanics, Vol. 1 (North-Holland, Amsterdam, 19’74). [16] J.B. Coon, N.W. Naugle and R-D. McKenzie, J. Mol. Spectry. 20 (1966) 106. (171 S. Bell. J. Phys. B2 (1969) 1001. [18] R.M. Wilcox, J. Math. Phys. 8 (1967) 962. [19] T. Ueda and T. Shimanouchi, J. Chem. Phys. 49 (1967) 470. [20] P.D. Mallinson and A.C. Robiette, J. hlol. Spectry. 52 (1974) 413. [X] R.A. CresweU and 1-M. Mills, 1. Mol. Spectry. 52 (1974) 392. [22] S. Crafn, V. Crecci and G. Turchelti, Nuovo Cimento 4B (19.71) 313.
[23] A. ErdClyi, W. Magnus, F. Oberhettinger and F.G. Tricomi, Tables of integral transforms (McGraw-Hill, New York, 1954). [24] I.S. Gradshteyn and I.M. Ryzhik, Tables of integrals, series and products (Academic Press, New York, 1965). 1251 R.K. Wehner and D. Baeriswyl, Physica 81A (1975) 129. [26] S.I. Ch+ and D. Sielman, J. Chem. Phys. 39 (1963) 545. [27] S. Bell, R. Davidson and P.A. Warsop, J. Phys. B 3 (1970) 113. [28] F.T. Hioe and E.W. Montroll, J. Math. Phys. 16 (1975) 194.5. [29] R.H.C. Hclleman and E.W. Montroll, Physica 74 (1974) 22.
.
:
. . .-‘.‘ UPPER AND LOWER ,ODNDS TO PARTITION FUNCTIONS -..: FGR MOLECULES : WITH-DOUBLE MINIMUM POTENTIALS .i._ __ .- .I. GRONHOLZ and W. WITSCHEL Departmht df C’henii.r&, UniirersiIy of Ulm. Germany .: Rx&d
i2 September 1977
Up&r and lower bounds to the viirational partition function 4 for molecules with double minimum potentials are derived. Both the exact lower (Gibbs-Bogoliubov) and upper (Golden-Thompsck) bound to q are evaluated analytically for the harmonic oscillator perturbed symmetrically or asymmetrically by a gaussian barrier. For the quadratic-quartic oscillator only the lower bound can be evaluated analytically, whereas the upper bound leads to a strongly divergent series resisting conventional summability techniques. In this case the classical partition function is used as an upper bound.
I. Introduction In a recent review on statistical methods for calculating thermodynamic functions, examples were given for molecules with double minimum vibrational potentials .[l] . Similar potentials are also interesting for hydrogen bonding and proton transfer [2,3] and for structural phase transitions in solids [4]. The customary technique is the separation of the remaining vibrational degrees of freedom and the numerical solution of the one- or two-dimensional Schriidinger equation either by direct solution
of the differential
tion [S] or by matrix dlagonallzation.
equa-
The molecular partition function q can then be evaluated by summation of the Boltzmann-factors calculated from the numerical eigenvalues. From a theoretical point of view it is more satisfactory to get analytical results which can be discussed generally. For a conventional force field the partition function can be obtained by means of thermodynamic perturbation theory leading in each order to simple results in ternis of exponential and hyperbolic functions [6]. As perturbations which cause the double minimum character of the potentialcannot be considered as small, different non-perturbatlonal.meth-odsmust beused. -Variational methods are non-p&turb&ive and have found wide.application in statistical mechanics. The recent reviews, by.Girardeau and Mazo [7] and
by Falk [8] summarize the mathematical foundation and applications of upper and lower bounds. Feynman [9] applied the lower bound to a linear oscillator with cubic anharmonicity which seems to be the only application to molecular vibrational problems. In this article we derive upper and lower bounds to the partition function q for the double minimum potentials mentioned. We shall not discuss again the mathematical background especially of the GoldenThompson inequality but refer to recent publications on this topic. The potentials behave well so that no principal difficulties arise. Therefore, only a brief sketch of the inequalities and the formalism to be used will be given in section 2. Section 3 contains the derivation of the bounds. The results are.evaluated numerically and discussed in section 4. The paper closes with a conclusion in section 5.
2. Variational principles and computational
methods
As the derivation of the following principles can be found in the literature [7-l I], we omit a repetition and just state the Gibbs-Bogoliubov and the Golden-Thompson inequality for the lower and upper bound to the partition function q. Consider a hamiltonian H, iiihich can be separated into
-’: .. _-: :
where fit, is not necessarily small. The Gibbs-Bogoliubov inequality represents.an upper.bound for the Helmholtz free energy A A=-~-llnq&40+(H1+r,
(2)
where 0-l = kT. A0 = -/1-l In-q0 is the_HehnhoItz free energy corresponding to HO, and (HI jn means - Tr{fil exp (-pfio)] _ @I& = (3) Tdexp t--P&>> Taking exponentials on both sides of eq. (2) yields the fohowing inequality for the partition functions 4 and qo q ~Tr{exp(-/I@)}
>q,-, exp (-@$~,)
=qL,_
(4)
which is the starting point for the calculations of lower bounds to ~7(the index “L” signifies “lower bound”). The Golden-Thompson inequality is 4 G qo(exp (-PHI PO = qu,
&xp[i(&
+ yj)] >n = exp (-$o~~ci>exp i-+y2y),
(10)
from which we .ob tain (11)
with (5)
(6)
In addition to the reviews mentioned above, we refer to the work of Breitenecker and Griimm [12] and of Ruskai [13] , dealing with the mathematical foundation of the Golden-Thompson inequality. Combining the inequalities, one has 4L =Q=Qu.
so that *r contains tberemainmg potential-energy terms.- This choicehas. the $dkta&; :&tit the &IcuIation of the &eragesm.@r>Oand
(exp (cZ))~ = exp <$fx2r-(>.
which will be used to compute upper bounds to 4. Taking Iogarithms on both sides one has a lower bound to the Heimholtz free energy A A >A, - fl-‘(exp (-&YL)J,.
’
(7)
The separation oftiinto two parts is fairly arbitrary, but from a practical point of view,.the number of possibilities is limited by the methods for evaluating. the traces in a suitable basis. If, for_example, one takes the kinetic energy operator for Ho and the potential energy operator for HII the averages can be calculated in the basis of plane waves with appropriate boundary conditions. Golden [lo] identified this qu as the cIassicaI partition function, qd- Therefore, one has another useful inequality
which wih be used to calculate 4” for the quartic OScihator. In ah other cases treated here, the appropriate
p = coth($&).
(12)
In the following we abbreviate /9zO by e. Formula (11) is the master formula for the subsequent calculations. It will be applied to the averages for the special cases we are interested in.
3. Calculation of the thermal averages 3.1. Potentials with gaussian barriers T&o types will be treated: ti‘asymmetric minimum potential with hamiltonian (I)
double
fi=fi0+&X2exp[-y(Z-a)2],
(I31
finding application in the theory of the hydrogen bond [2,3], and a symmetric double minbum potential with hamiltonian. (11) fl=Ho
+.(Yh~exp(-#),
.’ . . .
(14).
which is standard in the. molecular speckoscopi- of 1 large amplitude motions 116,171 .-Prelimbiaj c&u-lations using (II).as a model system in structtirz&ph&e transitions are prom&g. The principal forni of &se potentials is given in figs::l:and 2 bikmaji-v&y cotisiderably depending on the ~dten&d parameters.’ .. The thermal averages corres@%&ng to.(I) are given
287
.I. Gronholz, W. Witschel/Upper and lotier bounds to the molecularpartition function
power ofi appears in the exponential, and our master formula (11) cannot be applied directly. This difliculty is solved using the F,!urier-transformation exp[-r&i - A)2] = (4nn y)-lj2 +Ol X exp (--IrA-?/4ny) exp (i&) d7, s --m
(17)
which is a special applickion of Weyl’s prescription for quantizing classical functions, discussed extensively by Wilcox [18]. Formula (17) is introduced in eq. (16). The order of integration and trace formation is interchanged and Bloch’s theorem is applied -2 -1 COORDINATE
i
(clim&ionl&l
Fig. 1. Potential enera V(X) of hamiltonian (I). V(X) is the .sum of the harmonic oscillator potential V,(X) = (r!n/2)@ and the gaussian barrier V,(X) = fi 0 CLexp [ -T(X-A)~ ]_ Potential parametersare: (L= 2.7 = 1, A = 0.5.
(exp [-KY@? - A)‘] lo = (4 71~y)-l12 +m X exp (-h-A--?/4rry)(exp (iTJ@, dr s -m = (4 7111 y)-lj2
by
= (1 +izyp)- “’ :;
Gljo = k~iX’2exp [-r(_%! -
A)2] j,,
(15)
and (exp (-@,)+,
7
exp (-irA-T2/4rzy
.
- T2p/4) dr
[-A” /zy/( 1+izyr_1)]_
(18)
Combining this formula with eqs. (4) and (5); the upper and lower bounds to 4 for case I are
.
= (exp {-m exp [-~(2
4:’= q. expC--ad1+w)-~‘~exp[-A2d( 1+w)l 1, (1%
- A)2] }jo
=n$oy(exp[-rry(2-A)2])o,
(16)
*$I
=qo
5 I29
n-0
(1
t.yfi)-1’2
n!
X exp [-A2n-yl( 1 + JZYP)].
so that we are left with the calculation of (exp[-ny(* - A)2] jO_Unfortunately, the second
(20)
The corresponding formulae for the symmetric case II are obtained as special cases by putting A = 0. Thus, we get 4y
= q.
($0 = q.
ewI-4+rd-1121,
(21)
5 k$n.
(72)
( ltay/L)-“2.
n-0
The series for the upper bounds are convergent for any e, OLand 7 (7 > 0). The situation is less satisfactory for the quartic osci!lator as will be seen in section 3.2. COORDl&fATE
X (dimensionless1
Fig. 2. same as fig. 1, except that LI = 0 [hamiltonian (II)] in thgcase.
3.2. Quadratic-quartic
oscillator
The quadratic-quartic oscillator is of wide interest in quantum field theory, solid state physics, and statis-
288
J. Gronf~olz, W. .Wiisc?ze!lUpperand rower bounds to the moIecu~a~~artiti&fanHim
.tical mechanics. In molecular physics it is used to describe the strongly anharmonic puckering vibration of small ring molecules like cyclobutane [ 19] and oxetane [20,21] _Its hamiltonian is H=ti0m2(f&f4~).
(23
If& = -f , we have a pure quartic potential (case III) and iff? = 0, we simply add an z4-term to the harmanic oscillator potential (case IV). The averages
(24)
ff4G?>,),
and (exp (&I)),
= (exp[--ECf2J?2 tf,J?)]
lo
:
-~
.-.
Instead of expanding exp (-&I), one. may try to transform exp (-@gl) into anotherexpression, which can tie.averaged by means-of ouimasier formula; Using the Weyl prescription,is.given by Wilcox .[18], ._ -. we have to cair$late the Fourier integral g(~)
=$. J
exp I-eV,X2
’
+&x4) exp (-irx)dx.
-a.
: (30) However, the Fourier transform of the function exp (-X4) could not be found in the tables of Fou:. rier transforms [23] and, theiefore, this m&hod fails, too. From the above results it is clear, that we cannot average the operator function exp[-eCf2p2 +f,J!‘!)] using Bloch’s theorem. As the classical partition function qcl was shown to be an upper bound to q, we
calculate 1 %l=G (W
The expressions dm>u are deduced from our master formula, if we‘expand the exponentials and compare equal powers of (Yon both sides of (11). We get
(26)
i-a
ss
exp (-fed)
--m
X expf_-_4e[(1+2f2)X2 +2f,X4]}dPdX,
~cl=(~7r/e)1’2[(i X e’ K1,4(4
Inserting this into eqs. (4) and (5), we have qL = q. exp ]--E($ f2c1 f %f,&l
,
(27)
and
(31)
which can be separated into a translational and a configurational part. Both integrals are available from an integral table [24] . The result is +2f2)/2f4]112 -
Q4W1,
(32)
where z is equalto ~(1 + 2f2)*/(32f4). The 1, are modified Bessel functions of first kind and fractional order V. In the case of a pure quartic oscillator, the above expression simplifies to &II) = aI12 2-r f,-‘/4&~(3/4)-~_
x
(4n - 2/c)!
2n-k
(2n - k)! 42n-k ’
-
(28)
(33)
Eq. (32) is not new. Using confluent hypergeometric functions instead of the I,, it was applied by Wehner and Bear-My1 [25] to calculate the zero field susceptibility of a system of classical anharmonic oscillators.
For case IV, qr, reduces to
(4’*Y
3.3. Improvement.of
(2%
>z!(~,z)!42n Unfortunately, the inftite series in eqs. (28) and (29) diverge-for any E, f2 and f4 (E,f4 # 0). Customary summation methods, which were successful for the Rayleigh-SchrEdinger perturbation series of the quartic oscillator [22] , fail as the series are too strongly divergent.
the botinds
I Our bounds can be improved by the.introduction of variational parameters using,a suggestion of Feyn-
man [9]. . The hamiltonian is separated into fir, and tiI_ &=r-I, +fil_
(34)
A potential energy term TA =h a’(.? - A’)2 is added
J. Gronhols~W. Witschel/Upper and lower bounds to the molecular partition function
and subtracted (35) and-the aver&cs~(fjl - F$) and (exp [-/3(kl - &)]) are performed, using the density.operator corresponding to (fro f fh)_ The resulting expressions can be
Table 1 C&es I and IL Dependence of the percent difference between upper and lower bound, 070, on the potential parameters (Y, ‘y,and A. D% is related to the arithmetic mean $7” + qL) and given by eq. (37). fin/W = 1
varied under the conditions 3qll(cr)/&
a0
(a = !2’or A’; 7 = U or L.).
(36)
For symmetric potentials one finds, as expected, Abptimum = 0. In all other cases considered here, the extremum conditions (36) lead to transcendental equations which cannot be solved analytically, and
therefore, the variation has to be performed numerically. This is ti disadvantage, as we are mainly interested in closed results. Numerical work is postponed to ;1 later paper, dealing with coupled osciilators. Only the qtiartic oscillator will be given as an example, because in this case the improvements from the variational procedure are most significant_ The expression to be optimized with respect to CL’is then given by (27) if we putf2 = -l/2 and C2= .Q’.
4. Numerical
results
4. I. Po terrtials con taikng gatcssian bakers
In tables 1 and 2 the difference between upper and lower bound, D = q u - qL, is used as a measure for the accuracy with which we can estimate 4. combining qu and qL. Related to the arithmetic mean f(qu+qL), we express D in percent D% = 2OOD/(q,
+ qL).
289
2 2
2 6
0 0
2 6
10 2 6 10 2 6 10 2 2 2 10
0 0 0 0 0 0 0 2 4 2 4
-6 6 10 10 10 2 2 10 10
11.88 9.31 7.79 63.42 47.58 39.59 90.86
0.4034 0.5621 0.6302 0.0713 0.1929 0.2719 0.0126 0.0667 0.1173 0.6572 0.9298 0.1447 0.3506
76.58 66.91 7.12 0.62 62.31 38.12
0.5122 0.6774 0.7368
0.3186 0.5430 0.6283 0.2632 0.4988 0.5917 0.7580 0.9413 0.6231 0.7827
bounds are best for small e-values, i.e. high tempera-
tures and/or low vibrational frequencies. Finally, in table 3,qU and qL are compared to a partition function 4 calculated by summation of Boltzmann-facts,rs with numerical eigenvaules from the literature. The comparison shows that qu is always cIoser to CJthan qL. This means that the upper bound is a better approximation to 9 than the lower bound. 4.2. Quark-quadratic
potettriol
Though our results apply for any combination
of
f2and f4 provided f4> 0, we confine our numerical
(37) Table 2.
Table 1 presents LEG and the upper and lower bound to the partition function 4 for cases I and II. E is kept fEed (E = 1) and &potential parameters (Y,7 and A
CaseI.
are varied. One sees, that II% decreases with decreasing CYand increasing y and A. Small a; large y and A stand for a “small” I?l in the sense of ordinary perturbation theory. This means, our bounds improve with decreasing magnitude of the perturbation. The next question is, wh>ther qu and qL approx@ate q better for high or for low temperatures. The
paramctcrs:
answer can be seen from table 2, where the relative difference of the bounds is listed again for various evalues and a given perturbation. It appears that our
Dependence of the percent difference D% between upper and iower bound on E = ML//CT. D% is related to the arithmetic mean &J + 4~) and given by eq. (37). Potential OL= Y = A = 2
(0
E
QL
D%
Q{f)
0.2
4.646 2.107 1.281
0.4 1.5 3.2 5.1 7.1 16.2 33.0
4.680 2.170 1.364 0.980 0.758 0.329 0.104
0.4 0.6
0.8
0.884
1.0 2.0 4.0
0.657 0.237 0.053
:.
290
J: Gronholz, W. iVitschei&pper and Iower bounds to t&molecular
: .’ : -.. .’ f :. I_:. .. -:_.. .:
p&titio~ function-,
..
Table 3
-.
L- ,.
--‘_
-_:
COIIIPirkm of the upper and lower bound to the partition function q for cases I and IL q is calculated by summmg~u~ the
(I, 11) qL
A
Y
0.289
19.05
0
0.578 0.868 5.844
9.53 6.35
0 0
29.440 29.440
0.025 0.025
29.440
0.025
0.159
0.9177 0.8473 0.7655 0.6211 0.4123 0.3916 0.3809
0 -0.505 -0.429 -0.382
x x x x
:.
1c+ lo-l2 1O-‘2 lo-‘2
.. ,(r,IIj
~.
0.9190 0.8548 0.7845 0.1137 0.8136
a) a) a) I x 10-t x lo-‘t .0.6235.x lo-” 0.5333 x IO-”
:.- (I,II) 4rJ ._ :
... :
.‘-:. 0.9209. 0.8614
:. -:
.I.
. b) b) b) b1
0.7969. : 0.1266 x 10-t 0.1440 x lo-to
.-0.1095 x io-‘0 ‘0.9469 x lo-”
a) Eigenvalueswere taken from Coon et al. [16]. b) Eigenvalucs were taken from Flanigan and de la Vega [2] -
investigation to the quartic oscillator because its eigenvalues are known best [26,27] and the computation of the modified Bessel functions [l/4 in (32) is superfluous. In table 4 we list qrD and the corresponding “variationaIly improved bound”, called qzft3. The upper bound is given by the classical partition function qcl. qcl approaches 4 with decreasing E. This is exactly the high-temperature behavlour to be expected for a classical partition function_ In contrast to qcl, the lower bounds approximate 4 best in a medium e-range for E = 1. The variations
from case I, II and also IV (though not treated here) and arises from the fact that the hamiltonian (III) does not contain any harmonic contribution to the potential energy. The best result for &tl deviates approximately 3% from the numerical 4.. 4.3. Summary We summarize the results. In cases I and II our bounds are quite accurate, if the height of the gaussian barrier is smaller than or comparable to the.Ievelspacing of the corresponding harmonic oscillator (i.e. (Y< 1) and if E < 1. For larger perturbations and evalues the results are still useful but less accurate_ For a purely anharmonic potential as m case III, we cannot abandon the numerical variation, because the lower bound without variation gives very poor. results. On the other hand, it is advantageous, that qLv gives good estimates of 4 even. for extremely strong anharmonicity by,a minimum of computational effort. The numerical variation is very simple and can be performed-even on8 programmable pocket calculator. The upper bound, q&s only useful for. high temperatures, i.e. for E Q 1. -_
lead to considerable improvements, especially for high temperatures, where qLv and qL differ by a factor
of IO*-103. The lower bound without variation is nearly useless in this case. This situation is different Table 4 Comparison of upper and lower bounds to the partition function of the quartic oscillator, q(I1ll. &ll) and #1 arc givenby (27) and (33). E=trR/kT. Forsfff), see text l
(III) qL
0.2
0.004
0.4 0.6 0.8 1.0
0.877 0x939 0.2579
2.0 4.0 6.0 8.0
0.2847 0.2251 0.0774 0.0232 0.0067
(III)
q(III) al
,(fU c-_
2.7182 1.5892 1.1469 0.9002 0.7386 0.3580 0.1153 0.0390 0.0132
2.8450 1.6593 1.1943 0.9350 0.7654 0.3693 0.1204 0.0416 0.0144
2.8757 1.7099 1.2616 1.0167 0.8600 0.5113 0.3041 0.2243 0.1808
4LV
5. Conclusion--
:
:
1.
:: We applied variational techniques to. the partition: functiodqf one-diniensionid strongly-anh&mon& model systems, finding application-iii molecular physics, but also m-other.fieIds__: ‘: 1. :-.~’.I.,. -: .I- .- ‘.. : . The present work’airnkd $t two pomts: .First,: to :. .~
a) Calculated from eigenvahtes of Bell et al. [27] _
I
.:__
. .
i
._
-,_
_:
--..
_
_I.G~o&olz. W. WitscheljUpper and iower bounds to the krolecutar partition finction
formulae ‘for the one-dimen&I& double minirnUni problems, second, to test the appl&bi& of variational methods t6. the partition function of coupled tihamlonic 6scillators. The calculation of eigetivalues and eigenfunctions [28] and the general b&aviour [29] of such systems is a field of active reSearch with many open mathematical questions_ For the solution of chemical problems, especially in the field of reaction kineitcs, it is necessary to find good approximations to the partition function. _The satisfying results for one-dimensional strongly anharmonic systems encourage calculations for coupled oscillators.
‘derive u&Xul anal&l
Acknowledgement Constructive remarks on an earlier version of this paper by Professor R.A. Sack are gratefully acknowledged. We also wish to thank the “Fonds der Chemischen Industrie” for financial support.
References [I ] S.G. Frankiss and J.H.S. Green, in: A Specialist Period[21 131 i4j [S] [6]
ical Report, Chemical Thermodynamics, Vol. 1 (The Chemical Society, London, 1973). M.C. Flanigan and J.R. de la Vega, J. Chem. Phys. 61 (1974) 1882. R.W. Redding, J. Mol. Spectry. 62 (1976) 8. H. Beck, J. Phys. C9 (1467) 33. J.W. Cooley, Math. Comput. 15 (1961) 363. W. Witschel, 6. Grosswendt and J. Bohmann, PTBMitteihmgen 87 (1972) 3.
291
[7] M.D. Girardcau and R.M. Mazo, Advan. Chcm. Phys. 24 (1973) 187. [8] H. Faik, Am. J. Phys. 38 (1970) 858. [9] R.P. Feynman, Statistical mechanics (Benjamin, New York, 1972). [lo] S. Golden, Phys. Rev. B 137 (1965) 1127. [ll] C.J. Thompson, J. Math. Phys. 6 (1965) 1812. [12] M. Breitcnecker and H.R. Griimm, Commun. Math. Phys. 26 (1972) 276. 1131 M.B. Ruskai, Commun. Math. Phys. 26 (1972) 280.
[14] F. Bloch, Z. Phgsik 74 (1932) 295. [lS] A. Messiah, Quantum mechanics, Vol. 1 (North-Holland, Amsterdam, 19’74). [16] J.B. Coon, N.W. Naugle and R-D. McKenzie, J. Mol. Spectry. 20 (1966) 106. (171 S. Bell. J. Phys. B2 (1969) 1001. [18] R.M. Wilcox, J. Math. Phys. 8 (1967) 962. [19] T. Ueda and T. Shimanouchi, J. Chem. Phys. 49 (1967) 470. [20] P.D. Mallinson and A.C. Robiette, J. hlol. Spectry. 52 (1974) 413. [X] R.A. CresweU and 1-M. Mills, 1. Mol. Spectry. 52 (1974) 392. [22] S. Crafn, V. Crecci and G. Turchelti, Nuovo Cimento 4B (19.71) 313.
[23] A. ErdClyi, W. Magnus, F. Oberhettinger and F.G. Tricomi, Tables of integral transforms (McGraw-Hill, New York, 1954). [24] I.S. Gradshteyn and I.M. Ryzhik, Tables of integrals, series and products (Academic Press, New York, 1965). 1251 R.K. Wehner and D. Baeriswyl, Physica 81A (1975) 129. [26] S.I. Ch+ and D. Sielman, J. Chem. Phys. 39 (1963) 545. [27] S. Bell, R. Davidson and P.A. Warsop, J. Phys. B 3 (1970) 113. [28] F.T. Hioe and E.W. Montroll, J. Math. Phys. 16 (1975) 194.5. [29] R.H.C. Hclleman and E.W. Montroll, Physica 74 (1974) 22.