Alternative perturbation method for the molecular vibration–rotation problem

Alternative Perturbation Method for theMolecular Vibration–Rotation Problem

P. CASSAM-CHENAI,1 J. LIEVIN2

1Laboratoire d’Etude Theorique des Milieux Extremes, CNRS-OCA-UNSA, Physique-Recherche,Faculte des Sciences, Parc Valrose, F-06108 Nice, Cedex 2, France2Universite Libre de Bruxelles, Laboratoire de Chimie Physique Moleculaire, CP 160/09, 50 Av. F.D.Roosevelt, B-1050 Bruxelles, Belgium

Received 25 May 2002; accepted 8 January 2003

DOI 10.1002/qua.10556

ABSTRACT: This article introduces an alternative perturbation scheme to findapproximate solutions of the spectral problem for the rotation–vibration molecularHamiltonian. The method is implemented for the Watson Hamiltonian and applied tomethane. The complete J � 0 spectrum of this penta-atomic molecule of atmosphericinterest is calculated up to 9200 cm�1 in a purely ab initio fashion. Then, the rotationalspectra of the vibrational ground state is obtained up to J � 18. The largest relativeerror is 2.10�5 (for the highest J � 18 level) after scaling of a single parameter. Withoutscaling the accuracy on the rotational levels is limited by that of the ab initio equilibriumbond distance. The convergence of our results is analyzed with respect to the differentparameters involved in our approach. The important concept of vibrational mean-fieldconfiguration interaction is introduced. © 2003 Wiley Periodicals, Inc. Int J Quantum Chem93: 245–264, 2003

Key words: rovibrational spectra; perturbation theory; ab initio calculation; methanevibrational levels; methane rotational levels

1. Introduction

I t is remarkable that the diversity of numericalapproaches available to the quantum chemist

through widespread computer packages to calcu-late ab initio electronic or vibrational or rotationalwave functions decreases like the order of magni-

tude of the energy of the different spectral domains.As a matter of fact, the current distributions of themajor quantum chemistry codes usually includemany numerical methods for electronic calculationsbut do not go beyond the crude harmonic vibratorrigid rotor model. In recent years, the efforts tosolve the Schrodinger equation for the vibrationalmotion of polyatomic molecules have intensified,often triggered by astrophysics and atmosphericsciences. However, general computer codes im-proving on the harmonic oscillator approximation

Correspondence to: P. Cassam-Chenaı; e-mail: [email protected]

International Journal of Quantum Chemistry, Vol 93, 245–264 (2003)© 2003 Wiley Periodicals, Inc.

and able to treat a polyatomic molecule of an arbi-trary type are still scarce and far from enjoying thepopularity of the electronic-oriented packages. Themost cited are arguably SURVIB [1], SPECTRO [2],NIVELON [3], ANHAR [4], and MULTIMODE [5].Even fewer deal with the rotational motion of themolecules.

The purpose of this article is to introduce a newgeneral method to solve the Schrodinger equationfor the rovibrational motion of polyatomic mole-cules. A brief review of the other purely ab initiomethods that have been implemented as generalcomputer codes is provided to highlight the noveltyof our approach. We will not address such methodsas the empirical calibration of rotational constants[6] nor the averaging of the �-matrix elements onvibrational wave functions or by means of quantumMonte-Carlo calculations [7] because they do notpropose a general approximation scheme converg-ing to the exact solution. We will not describe eitherthe methods that are not compatible with the use ofthe Watson Hamiltonian because, so far, they haveonly been amenable to system-specific implemen-tations.

It must be noted that the Watson Hamiltonianitself is not completely general and that two dis-tinct expressions must be employed according tothe molecule being linear [8] or nonlinear [9], butthis can be easily accommodated for in a singleprogram. Real problems arise when, for instance,a floppy motion distorts a nonlinear molecule toa linear geometry. Then, a singularity occurs inthe Watson Hamiltonian and the use of a modi-fied treatment to extract the floppy vibrationalmodes [10, 11] becomes unavoidable. Even innonpathologic cases, judiciously chosen non-Car-tesian coordinates can be advantageous to ex-pand the potential energy surface (PES) but at theexpense of the loss of generality mentioned aboveand the introduction of kinetic operators whoseexpressions become rapidly intractable as thenumber of atoms grows [12]. Although some suc-cessful attempts at dealing with relatively largesystems have been reported, the larger systemsbeing often treated by using a clever interplay ofbasis set representations for the different opera-tors of the Hamiltonian, there is currently nogeneral code available for molecules with fivenuclei or more, as far as we are aware. A reviewon this topics can be found in [13]; more recentreferences are, for instance, [14, 15].

The method to solve the rotation–vibration prob-lem most easily applied to a general system is ar-

guably the perturbation approach based on an iter-ative use of contact transformations [16]. Thistraditional perturbation method first formulated byNielsen [17] for a general polyatomic molecule(with no three-fold degenerate vibration) has beenreviewed in a clear fashion by Mills [18] or Alievand Watson [19] and will only be briefly presentedin Section 3.1. In this approach, the zero-orderHamiltonian is taken to be the harmonic oscillatorvibrational Hamiltonian (plus sometimes the rigidrotor Hamiltonian; see, e.g., Ref. [20]), the rest of theHamiltonian being a perturbation whose dominantterm is of order (a/l ), where a is a typical vibra-tional amplitude and l a typical bond length. As aresult, large amplitude vibrations for H-containingmolecules can cause the perturbation series to con-verge slowly [21].

A less conventional perturbation approach con-sists of taking in the zero-order Hamiltonian all theterms free of rotational operator [22–24]. The per-turbation is then of order (�rot/�vib) where �rot is atypical rotational frequency and �vib a typical vibra-tional frequency. This is usually an order of mag-nitude smaller than (a/l ); thus, one can expect theperturbation series to converge faster than that ofthe traditional approach. It is this perturbationscheme based on a somewhat unusual partitioningthat we investigated in this article. The zero-orderHamiltonian, unlike the harmonic vibrator Hamil-tonian, is not trivially diagonalized. It has beentreated variationally with a scheme of symmetry-adapted contractions reminiscent of the one used in[25] but totally flexible in this work. This new vari-ational method, called vibrational mean-field con-figuration interaction (VMFCI), proved particularlyeffective.

The rest of the article is organized as follows:In Section 2 we introduce the main variationaltechniques that have been employed for the rota-tion–vibration problem of a general molecule. Formolecules containing more than four nuclei, apure variational treatment becomes a hard task,in particular for J large, and one is led to elimi-nate some elements in the Hamiltonian matrix [4,26, 27] or use a perturbation scheme. Two differ-ent perturbation methods applicable to a generalmolecule are presented in Section 3: the tradi-tional Van Vleck method and a newly imple-mented method. In Section 4 we apply the newmethod to CH4 and analyze the convergence ofour results.

CASSAM-CHENAI AND LIEVIN

246 VOL. 93, NO. 3

2. Variational Resolution of theWatson Hamiltonian EigenvalueProblem

The general form of the rotation–vibration Ham-iltonian for an arbitrary choice of coordinate systemin the frame of the Born–Oppenheimer approxima-tion [28] has been given by Pickett [29]. It has beennoted that choosing a convention for the body-fixedframe (such as the Eckart convention) amounts to achoice of gauge [30, 31]. In this work we only fo-cused attention to the mass-weighted Cartesiannormal coordinate system in the Eckart frame. Ithas been demonstrated that, whether one startsfrom the “all motions” Cartesian Hamiltonian andextracts its rovibrational part [32] or, following Dar-ling and Dennison [33], applies the general quanti-zation rules of Podolsky [34] to the classic Wilsonand Howard [35] rovibrational Hamiltonian, oneobtains after simplifications the Watson Hamilto-nian [9].

In this section we give an overview of what weconsider the main, general variational methods thathave been applied to the Watson Hamiltonian sofar. The methods of this section (and also of Section3) do not specifically require the use of normalcoordinates. However, for the sake of simplicity, wewill restrict their presentation to a Hamiltonian ex-pressed in terms of vibrational normal coordinates(Qi)i�1,. . .,n and their conjugate momenta (Pi)i�1,. . .,n.The J � 0 part of the Watson Hamiltonian or moregenerally its representation in an appropriate rota-tional basis set assumes this form [36].

2.1. SELF-CONSISTENT FIELD METHODS

As for electronic calculations, it is usually inter-esting to begin with an “independent” vibrationalmode model, in which the n-mode wave function isexpressed in terms of single-mode vibrational func-tions (or modals), (�vi

(Qi))i�1,. . .,n, where vi is aninteger labeling the modal of the ith mode. A trulyindependent model where the force constants thatcouple several modes together are completely ne-glected can give poor results. The so-called uncou-pled anharmonic oscillators (UAO) approach [37]can even predict a wrong anharmonic behavior be-cause the intermode anharmonic constants canhave larger quantitative effects than intramode an-harmonic constants [38]. However, a vibrationalself-consistent field (VSCF) ansatz avoids thesesorts of discrepancies and provides a good starting

point in general for about the same computationalcost. The VSCF method [21, 39] proceeds as followsfor an n-oscillator system represented by a vibra-tional Hamiltonian of the general form

H � �i1

h1�Qi1, Pi1� � �i1,i2

h2�Qi1, Pi1, Qi2, Pi2�

� · · · � hn�Qi1, Pi1, Qi2, Pi2, · · · , Qin, Pin�, (1)

where the operator hp gathers all the terms involv-ing an explicit coupling between p oscillators, re-gardless of the actual origin of the coupling (kineticor potential energy). Note that an additional con-stant term h0 may occur but can be removed bysetting the origin of the energy unit at the zero-point energy (ZPE). The trial function used in theVSCF procedure is a single product of modals alsocalled a “vibrational configuration,”

�v1,v2,· · ·,vn � �i�1

n

�vi�Qi�, (2)

and characterized by the set of integers (v1, v2, . . . ,vn). This trial function leads to the system of VSCF-coupled equations, for j � 1, . . . , n,

�h1�Qj, Pj�

� ��i�j

�vi�Qi��H � h1�Qj, Pj�� i�j

�vi�Qi�� vj��vj�Qj� � 0, (3)

or, more compactly in self-explanatory notations,

�HjVSCF � �vj��vj�Qj� � 0. (4)

Each equation describes the vibrational motionalong a given coordinate in the mean field createdby the other oscillators. As in the electronic context[40], approximate solutions of the VSCF equationscan be obtained iteratively within the algebraic ap-proach by expanding each modal on a basis set ofappropriate analytic functions (�ri

(Qi))ri�1,. . .,Ni,

�vi�Qi� � �ri�1

Ni

avi,ri�ri�Qi�. (5)

MOLECULAR VIBRATION–ROTATION PROBLEM

INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY 247

The VSCF procedure then consists of iterativelydiagonalizing the matricial representations of the nVSCF operators Hj

VSCF of Eq. (4) in the Nj-dimen-sional vector space spanned by the (�rj

)rjbasis sets.

At convergence, the only single-mode functionsthat are optimized, in a mean field point of view,are those of the trial function, Eq. (2), i.e., charac-terized by the set of integers (v1, v2, . . . , vn). Oneobtains, however, in addition Ni � 1 virtual singlemode functions for each mode i. These virtualmodals are the building blocks of other n-modeconfiguration functions of the form given by Eq. (2)but for different sets of integers. They can be usedin further calculations, like vibrational configura-tion interaction (VCI), to introduce explicitly boththe inter- and intramode correlations.

The SCF approach has been extended in differentways to improve the accuracy of the VSCF solu-tions. Thompson and Truhlar [41] performed simul-taneous SCF optimizations of the modals and thecoordinate system. This idea was also exploited byother authors [42] (and references therein). Sch-wenke [43] extended the self-consistent procedureto a vibrational multiconfiguration trial functionusing a second-order algorithm. Culot and Lievin[44] also developed a vibrational multiconfigura-tion SCF method (VMCSCF) based on the use of thegeneralized Brillouin theorem [super-configurationinteraction (CI) algorithm] and showed [45] that thelatter approach was well adapted to implement thevibrational analog (VCASSCF) of the complete ac-tive-space SCF (CASSCF) method widely used inelectronic structure calculations [46, 47].

2.2. CONFIGURATION INTERACTIONMETHODS

The VCI method is the usual variational ap-proach for going beyond SCF. The trial function isexpressed as a linear combination of configurationfunctions �v1,v2, . . . ,vn

of the form given by Eq. (2):

�iVCI � �

�v1,v2,· · ·,vn�

c�v1,v2,· · ·,vn�,i�v1,v2,· · ·,vn. (6)

The modal basis set used to build the vibrationalconfigurations is arbitrary. However, a good choiceis usually the basis set resulting from a previousground-state VSCF calculation.

The matrix representation of the Hamiltonian,Eq. (1), in the vibrational configuration basis set(product basis set) must be diagonalized. The dif-ferent roots of this secular problem provide, in a

single calculation, approximate solutions for all thevibrational states subtended by the basis set. Thelabel i in Eq. (6) refers to the ith root of the matrix.Also, note that the H matrix can already be in ablock-diagonal form by construction, when the con-figuration functions are adapted to the symmetrypoint group of the molecule.

In parallel with the usual quantum chemistryterminology, a full-VCI calculation will involve allpossible configurations arising from the modal ba-sis sets. Full-VCI calculations are, however, limitedto small molecules, triatomics in practice, due to therapid blowup of the Hamiltonian matrix size withthe number of nuclei and the size of the modal basissets. Truncated VCI calculations are thus to be con-sidered, and an efficient solution adopted by differ-ent authors is the contracted CI or diagonalizationtruncation procedure [48, 49]. Separate Hamilto-nian matrices, built for limited sets of nuclear de-grees of freedom, are diagonalized and the result-ing set of eigenvectors is truncated on the basis ofenergy or individual mode quantum numbers cri-teria. A basis set truncation being carried out ateach intermediate step of a given contractionscheme, one can expect a drastic damping of thefinal VCI basis set size with respect to the full-VCIone. Different coupling schemes introducing pro-gressively in a hierarchical way the intermode cou-plings are, of course, possible.

2.3. MEAN-FIELD CONFIGURATIONINTERACTION METHOD

We developed in the CONVIV computer code[50] a flexible contracted VCI algorithm, whichtakes advantage of the mean-field approach de-scribed in Section 2.1 at each contraction step.

Consider an n-oscillator problem and a givenpartition of the n modes into q sets I1, I2, . . . , Iq, of,respectively, p1, p2, . . . , pq modes:

�I1, I2, · · · , Iq� � �i11, i2

1, · · · , ip11 ,

i12, i2

2, · · · , ip22 , · · · , i1

q, i2q, · · · , ipq

q �. (7)

For each contraction, Ij, we define a partial Hamil-tonian,

Hj � �i1�Ij

h1�Qi1, Pi1� � �i1,i2�Ij

h2�Qi1, Pi1, Qi2, Pi2�

� · · · � hpj�Qi1j , Pi1

j , Qi2j , Pi2

j , · · · , Qipj

j , Pipj

j �, (8)


248 VOL. 93, NO. 3

and the mean-field Hamiltonian

�Hj � � �Ik�Ij

�Vk�Qi1k, · · · , Qipk

k ��H � Hj�

�Ik�Ij

�Vk�Qi1k, · · · , Qipk

k �� Vj��Vj � 0, (9)

which is diagonalized on the full product basis set��Vj�Qi1

j , . . . , Qipj

j ��Vj with

�Vj�Qi1j , · · · , Qipj

j � � �i�Ij

�vi�Qi�, (10)

the superlabel Vj standing for (v1, . . . , vpj). Iterating

this scheme until convergence would lead to SCF-contracted solutions for the Vj modes. We havefound, however, that in practice two or three itera-tions before applying a truncation are largely suf-ficient in most cases. So, we term our method VM-FCI rather than VSCFCI because SCF convergenceis not systematically sought after.

An interesting feature of the contraction proce-dure is to properly deal with degenerate vibrations.It suffices to contract together all the degeneratecomponents to ensure the preservation of symme-try. This is in contrast with the traditional VSCFscheme, where the degenerate components are op-timized on their own, which usually results in anartifactual symmetry breaking.

Note that the VSCF of Section 2.1 and the VCI ofSection 2.2 are both particular case of VMFCI. Theformer corresponds to the partition (I1, I2, . . . , In) �({1}, {2}, . . . , {n}), the VMFCI step being iterateduntil a given convergence criterion is reached. Thelatter corresponds to the partition (I1) � ({1, 2, . . . ,n}) and is not strictly speaking “meanfield” becauseall modes are already included in the contraction.

When a given contraction step, (I1, I2, . . . , Iq), hasbeen performed, a new contraction scheme ( J1,J2, . . . , Js) can be designed provided that it corre-sponds to a coarser partition in the sense that s � qand

Ik � �I1, I2, · · · , Iq� Jl � � J1, J2, · · · , Js�

such that Ik � Jl. (11)

Trials and errors have shown us that, as in thetraditional VSCF, the best results at step x � 1 areusually achieved with the product basis set whosecomponents are obtained by means of a ground-

state VMFCI at step x, that is, by setting Vk � 0 inEq. (9). All the results presented in Section 4 corre-spond to ground-state VMFCI. The product basissets were truncated by setting an energy thresholdon the energy of one or several components or onthe sum of the energy of all the components.

The CONVIV program implements the VMFCImethod for a Hamiltonian of the form, Eq. (1),where all the hp are polynomial functions,

hp�Qi1, Pi1, · · · , Qip, Pip�

� �r1,· · ·,rp,s1,· · ·,sp

cr1,· · ·,rp,s1,· · ·,spQi1

r1Pi1

s1 · · · Qip

rpPip

sp , (12)

and with harmonic oscillator functions (solutions ofthe harmonic part of the potential) as modal basisfunctions, which is the more convenient choice todeal with a Hamiltonian in normal mode coordi-nates. The pure vibrational part of the WatsonHamiltonian can be cast in the form of Eq. (12) bylimiting the series expansions of the potential andof the � matrix to a given order (see Sections 3 and4).

3. Perturbative Resolution of theWatson Hamiltonian EigenvalueProblem

Perturbation theory has a long history in the fieldof the rotation–vibration problem. We referred to afew good reviews in Section 1. The aim of ourpresentation here is to highlight the differences be-tween the common form of perturbation theoryencountered in the field and the alternative methodthat we implemented.

3.1. VAN VLECK PERTURBATION METHOD

Van Vleck perturbation theory [16] consists ofthe repeated use of well-chosen unitary transforma-tions applied to an Hamiltonian matrix. The pur-pose of these transformations, called “contact trans-formations,” is to bring corrective contributionsdue to off-diagonal matrix elements of the originalmatrix on the diagonal of the new matrix. Theunitary character of the contact transformation en-sures the invariance of the predictions of quantummechanics in the new representation, hence theterms “canonical transformations” and “canonicalperturbation theory” also found in this context.



This approach has proved useful to both theore-ticians and experimentalists to relate the rotationalconstants of the effective spectroscopic Hamiltoni-ans to first-principle quantum mechanics constants.We will not address here the problem of retrievingforce fields and geometry parameters from the ro-tational constants and refer the interested reader tothe recent review by Herman et al. [13]. We willonly be concerned with the “direct” problem ofobtaining rovibrational spectra from force fields.

The first general formula based on the Darling–Dennison Hamiltonian were formulated by Nielsen[17] in the frame of Van Vleck perturbation theory.Later, Watson [9] simplified the Darling–DennisonHamiltonian and stimulated new investigations onthe relationship between the complete molecularHamiltonian and the effective rotational Hamilto-nian corresponding to a given vibration [18]. Theformula for the effective rotational constants de-rived from this perturbation theory were imple-mented in several computer programs [1, 2] andsystematic studies of the vibration–rotation interac-tion in polyatomic molecules became feasible [51].

Primas made the proposal to take advantage ofthe connections of this perturbative approach withLie algebra [52] and Sibert took up the challenge toapply Primas formalism to the field of the rovibra-tions of polyatomic molecules [53], followed byother authors.

Most of the implementations mentioned abovestart from a product of harmonic oscillator func-tions and rigid rotor functions as zero-order wavefunctions. However, a few exceptions [26, 54] usethe J � 0 vibrational Hamiltonian eigenfunctions inplace of the harmonic oscillator functions as in thealternative perturbation theory below. The majordifference with the latter approach lies in the par-titioning of the Watson Hamiltonian. The WatsonHamiltonian can be written as a sum [18],

H � H2,0 � H3,0 � H4,0 � H2,1 � H0,2

� H5,0 � H3,1 � H1,2 � H6,0 � H4,1 � H2,2, (13)

after truncation at sixth order of the PES and atsecond order of the � matrix,

� � �r�0

� �12

r

�r � 1�

� �k1,. . .,kr

Ie�1ak1Ie

�1 . . . akrIe�1Qk1 . . . Qkr, (14)

where Ie�1 is the inverse of the inertia tensor

I(Q1, . . . , Qn) at equilibrium and (ak)k the deriva-tives of the latter with respect to the normal coor-dinates,

ak � � �I�Qk

0

. (15)

The Hmn in Eq. (13) are the terms of the Hamiltonianof degree m in the vibrational operators and ofdegree n in the rotational operators. Their detailedexpressions are given in [18]. In the traditional ap-proach one follows Oka’s order of magnitude clas-sification [55], that is, that the perturbation expan-sion of the Hamiltonian with respect to the Born–Oppenheimer parameter � 1/10 assumes theform

H � H0 � H1 � 2H2 � 3H3 � 4H4 � · · · ,

(16)

where we have set

H0 � H2,0,

H1 � �1H3,0,

H2 � �2�H4,0 � H2,1 � H0,2�,

H3 � �3�H5,0 � H3,1 � H1,2�,

H4 � �4�H6,0 � H4,1 � H2,2�.

As we shall see, the perturbation expansion pre-sented in the next section corresponds to the expan-sion

H � H0 � 2H1, (17)

where this time

H0 � H2,0 � H3,0 � H4,0 � H5,0 � H6,0 � · · · ,

H1 � �2�H2,1 � H0,2 � H3,1

� H1,2 � H4,1 � H2,2 � · · ·).

3.2. GENERALIZED RAYLEIGH–SCHRODINGER PERTURBATION METHOD

We present here a generalized Rayleigh–Schrod-inger perturbative treatment for a Hamiltonian ofthe form


250 VOL. 93, NO. 3

H�X, Y� � H0�X� � �H1�X, Y�, (18)

that is, a Hamiltonian whose zero-order term hasno dependence upon certain coordinate and mo-mentum operators. This is, for example, the form ofthe Watson Hamiltonian [9] for nonlinear mole-cules if we set (in a.u.) X � [(Qi)i, (Pk)k] the vector ofthe vibrational mass-weighted Cartesian displace-ments from the equilibrium geometry and theirconjugate momenta and Y the vector of the Eulerangles in the Eckart frame and their conjugate mo-menta.

H0�X� �12 �

k

Pk2 � V �

12 �

��

�� 18 �

�

��,

(19)

�H1�X, Y� �12 �

��

�� 2��. (20)

In the equations above, V is the potential ofelectronic origin in the Born–Oppenheimer approx-imation, expressed as a function of the normal co-ordinates Qi, � is the 3 � 3 matrix whose seriesexpansion in terms of the normal coordinates isgiven in Eq. (14), � is the so-called vibrational an-gular momentum, which only depends upon theoperators in the set X, H0 is what is commonlyreferred to as the J � 0 Watson Hamiltonian, and �is the total angular momentum and the sole quan-tity depending upon the operators in the set Y.

There are two ways of considering a Hamilto-nian H as given in Eq. (18). First, one can say that Hdepends upon two sets of operators X and Y. Then,the standard Rayleigh–Schrodinger perturbationtheory can be applied, the nondependency of H0upon Y being just a particularity of the problem.Alternatively, one can take H as depending onlyupon the operators X and apply a generalized per-turbation theory, “generalized” because of the extradependency of H1 upon Y. As we are now going tosee, this second approach proves more fruitful thanthe first.

So, ignoring Y, we look for “x-eigenvectors”��(x) of H,

H�� x� � �� x�, (21)

of the form

�� x� � �0� x� � ��1� x� � �2�2� x� � . . . , (22)

associated with “x-eigenvalues” that are operatorsstill depending upon Y,

��Y� � �0�Y� � ��1�Y� � �2�2�Y� � . . . , (23)

that is, we only develop in power of � the part of thewave function depending on the set of spectralvalues, denoted by x, of the coordinate operators ofset X. In the case of the rovibrational problem, �� isa vibrational wave function and �� is the effectiverotational Hamiltonian associated with the vibra-tional level ��. This scheme resembles the adiabaticapproximation but, here, the “frozen” variables canbe the “slow” variables, e.g., the rotational variablesin the case of the Watson Hamiltonian. �� is nor-malized to unity,

�� 1, (24)

and its phase is chosen so that ��0�� be real (it isunderstood that the Dirac bracket means integra-tion over x only). In the usual method these condi-tions would be imposed to the complete wave func-tion of the system, which would also depend upony, the set of spectral values of the coordinate oper-ators of set Y.

At present we expand Eq. (21) up to a givenorder in � to obtain a reduced Hamiltonian ��(Y)depending only on operators of Y.

3.2.1. Zero Order

Equation (21) becomes

H0�0 � �0�0, (25)

which is just the usual eigenvalue equation forHamiltonian H0. Therefore, �0 is a constant inde-pendent of Y and if, in a usual Rayleigh–Schrod-inger approach, we were looking for solutions de-pending also on y we would obtain an infinitelydegenerate set of eigenvectors because their depen-dence upon y would not be determined by Eq. (25).

3.2.2. First Order

First, let us assume that �0 is nondegenerate. Atthis order, the reduced Y-Hamiltonian, �0 � ��1(Y),really depends on Y because �1(Y) is given by

�1�Y� � ��0�H1�X, Y��0�. (26)



It provides a new eigenvalue problem for the y-dependent part (i.e., the rotational wave function inthe case of the rovibrational problem) of the totalwave function:

��0 � ��1�Y�� y� � E�� y�. (27)

Each product �0(x)�(y) will be an approximate eig-envector of H associated with eigenvalue E.

Remark 1. In the usual approach, assuming thatthe solution has the factorized form, �(x)�(y), wewould have at zero order an infinite set (�0(x)�i(y))i

of degenerate functions as already noted. Thus,following the Rayleigh–Schrodinger perturbationtheory for the degenerate case and selecting somefinite basis set, at first order we would have todiagonalize the Hamiltonian matrix (��i�0�H1��0�j�)i,j.The scheme presented here allows for other possi-ble treatment of Eq. (27) than the total diagonaliza-tion.

Let us now consider the case of a degenerate set(or quasidegenerate with respect to the order ofmagnitude of the largest couplings) of orthonormalzero-order x-eigenvectors, (�0

�(x))1��. The pertur-bation may remove the degeneracy and a prioridifferent wave functions ��(y) are to be associatedwith each �0

�(x). However, we can obtain an eigen-value problem formally identical to Eq. (27) if �(y)instead of being a function of the Hilbert spacecorresponding to the y degrees of freedom, �y, isnow considered an element of the extended space,�y V ��, which is the tensor product of �y and a(spin-like) internal degree of freedom Hilbert space,�� (� is the field of complex numbers). In the caseof the rotation–vibrations, this amounts to consid-ering the rotational wave function associated with a�-fold degenerate vibrational level to be a functionwith � internal vibrational degrees of freedom.Then, �(y) can be decomposed into a direct sum of� components,

�� y� � �1��

�� y�, (28)

where ��(y) is a short notation for ��(y) V �, �referring to an element of the canonical basis set of��. A basis set of �y, (�i(y))i, induces in a naturalway a basis set of �y V ��, denoted (�i

�(y))i. Eachcomponent ��(y) can be expressed in the inducedbasis set as

�� y� � �i

ci��i

�� y�, (29)

and combining Eqs. (28) and (29) we obtain thefollowing decomposition of the wave function:

�� y� � ��

�i

ci��i

�� y�. (30)

The first-order corrective term of the reduced Ham-iltonian, ��1(Y), is decomposed into a double sum,

��1�Y� � � ��,�

�1�,��Y�, (31)

where

�1�,��Y� � ��0

��H1�X, Y��0�� t�,� (32)

is defined to be the tensor product of an operator��0

��H1��0�� acting only in the y-Hilbert space �y and

an internal degree of freedom operator t�,�

t�,��i� � ��,��i

�, @i, (33)

which changes the �-component of a wave functioninto an �-component with the same y-part (��,� isthe Kronecker symbol).

So, finally, the reduced Hamiltonian, [�0 ��1(Y)], will act on �(y) in the following way:

��0 � ��1�Y�� y�

� ��

�i

��0ci� � �ci

��0��H1�X, Y��0

��i�� y�. (34)

This Hamiltonian will partially lift the degeneracy.Consider, for example, the Watson Hamiltonian ofmethane that will be used in Section 4. The �4vibrational state of methane is three-fold degener-ate and carries the F2 irreducible representation ofTd symmetry. The combination of the J � 0 vibra-tional functions with the spherical harmonic func-tions of the Euler angles corresponding to J � 1(which carries the F1 irreducible representation)gives nine degenerate levels at zero order. The di-agonalization of the reduced first-order Hamilto-nian, Eq. (34), splits these rovibrational levels intofour states carrying the irreducible representationsA2, F2, F1, E in increasing order of wave number.

It is important to take into account this splittingat the higher orders of the perturbation theory aswell as the fact that the associated wave function


252 VOL. 93, NO. 3

can no longer be factored into an x-part and ay-part, i.e., a vibrational times a rotational wavefunction in the case of the rotation–vibration prob-lem.

3.2.3. Higher Orders

We begin again with the nondegenerate case.Let us denote by ( fk)k the other eigenvectors of H0(different from �0) associated with the set of eig-envalues (�k)k. We obtain an improved reducedY-Hamiltonian eigenvalue problem for �( y),

��0 � ��1�Y� � �2�2�Y�� y� � E�� y�, (35)

where

�2�Y� � �k

��0�H1�X, Y��fk��fk�H1�X, Y��0�

�0 � �k. (36)

We draw attention to the fact that a solution ofEq. (27) will not satisfy Eq. (35) in general, forthere is no reason for it to be an eigenvector of thesecond-order corrective term �2�2(Y). Equation(35) is a new eigenvalue problem in its own rightand, like Eq. (27), one is free to solve it by themethod of one’s choice. For instance, one cantreat the second-order term as a perturbationwith respect to the first-order equation, Eq. (27).However, it may be convenient sometimes to in-

corporate part of the operators of the second-order correction into the first-order Hamiltonian.Alternatively, one can perform a total diagonal-ization in a chosen basis set (�i( y))i of the wholereduced Hamiltonian.

In the case of the Watson Hamiltonian, the �2

operator polynomial expansion in terms of the an-gular moment operators has terms of degree 2, 3,and 4. After some rearrangements using the com-mutation relations of the angular moment opera-tors one can get rid of the terms of order 3 andmake appear the leading contribution to the quarticcentrifugal distortion constants and a correctivecontribution to the rotational constants.

In the degenerate case, because the first-orderperturbation theory would have given a wavefunction of the form ¥�,i di

��0�( x)�i

�( y), where thex-part and the y-part cannot be separated out, wehave to go back to the regular Rayleigh–Schrod-inger perturbation theory and apply this methodto the total problem. However, the partitioning ofthe Hamiltonian Eq. (18) still holds true and isexpected to lead to a rapidly converging energyexpression.

We stop here the description of the method,but it is clear that it can be extended to anyarbitrary order. We give below without furtherdemonstration the expansions of �3(Y) and�4(Y):

�3�Y� � �k,k�

��0�H1�X, Y��fk��fk�H1�X, Y��fk��fk��H1�X, Y��0�

��0 � �k��0 � �k��

� ��0�H1�X, Y��0� �k

��0�H1�X, Y��fk��fk�H1�X, Y��0�

��0 � �k�2 , (37)

�4�Y� � �k,k�,k�

��0�H1�X, Y��fk��fk�H1�X, Y��fk��fk��H1�X, Y��fk��fk��H1�X, Y��0�

��0 � �k��0 � �k��0 � �k�� 0�H1�X, Y��0�

� �k,k�

��0�H1�X, Y��fk��fk�H1�X, Y��fk��fk��H1�X, Y��0�

��0 � �k��0 � �k�� 1��0 � �k�

�1

��0 � �k��

� ��0�H1�X, Y��0��2� �k

��0�H1�X, Y��fk��2

��0 � �k�3 � �

k

��0�H1�X, Y��fk��2

��0 � �k��k

��0�H1�X, Y��fk��2

��0 � �k�2 . (38)

Equations (26), (36), (37), and (38) when applied tothe Watson Hamiltonian, Eqs. (19) and (20), givethose obtained by Makushkin in [23] (with somemisprints corrected) for his vibrationally averaged

rotational Hamiltonians. The general formula forthe nth-order perturbation corrections of Kato andthat of Bloch found in Ref. [56] can be used toconstruct the higher-order �n(Y).



Remark 2. A similar perturbative method is alsoapplied in a standard textbook to the diatomic case,with the summation over the vibrational states re-stricted to the states immediately above and belowthe zero-order vibrational state [57].

Remark 3. Another alternative to the traditionalscheme is Møller–Plesset perturbation theory. It isnot addressed here because it has only been imple-mented so far for the pure vibrational problem [58].The dimensional perturbation theory applied to therovibrational problem of triatomic molecules by Su-vernev and Coworkers [59] is worth mentioning aswell, but it has not been exploited to deal with thegeneral case, as far as we are aware.

4. Application to Methane

In this section we apply the perturbation methodof Section 3.2 to the main methane isotopomer aftera variational treatment of the vibrational degree offreedom, and study its convergence properties.

The choice of methane was dictated by physicalas well as technical considerations. Methane is aspecies of interest in astrophysics and atmosphericsciences. It has been detected in the interstellarmedium, the atmosphere of brown dwarfs, andmany objects of the solar system. These include inparticular Earth’s atmosphere, where methane is agreenhouse gaz, Titan’s atmosphere, where lifecould have evolved, Jupiter, where hot methanebands observed after the collision with comet Schu-maker–Levy 9 showed the limit of the spectroscopicdatabases [60], and comets themselves. The urge forimproving our knowledge of methane infrared (IR)spectroscopy has been strongly stressed by thebrown dwarf and extrasolar planet scientific com-munity [61, 62]. As a result, it has been the subjectof many studies both experimental and theoretical.However, it remains a challenging species becauseof the density of its spectrum (it is a penta-atomic;hence, it has 12 degrees of rotation–vibration) andthe numerous couplings (including importantFermi resonances and Coriolis couplings) that pre-vent a simple interpretation of its spectrum. Meth-ane provides also a good test for the convergence ofour perturbation method because it has relativelylarge-amplitude hydrogen motion. Many recenttheoretical studies are available for comparison [63–71]. Its high tetrahedral (Td) symmetry imposesmany degeneracies that are useful tools to diagnosebugs in computer programs. Finally, it is the sim-

plest spherical top molecule and all other types ofrotating systems (symmetrical top, asymmetricaltop) can be obtained by isotopic substitutions andstudied with the same PES.

Notations. Our vibrational calculations are donein two steps: First, the program WATCOR expandsthe � matrix at order n for H0 and l for H1, combinesthe Born–Oppenheimer potential of order m withthe so-called “Watson term,” and forms the J � 0Hamiltonian to be used by the program CONVIV.The latter solves the vibrational eigenvalue equa-tion using the VMFCI method and calculates theintegrals needed for the calculation of the H1 (andof observable like the dipole moment) matrix ele-ments. In the following the notation WmnlCpq willrefer to a calculation consisting of a WATCOR runcorresponding to the orders m, n, l mentionedabove, followed by a CONVIV run that includes allstretching modes up to the approximate quantumnumber p and all bending modes up to the approx-imate quantum number q. Further, we will set n tothe letter “x” when no rotational contribution isincluded in the vibrational Hamiltonian and l to theletter “y” to designate a WATCOR run where onlythe purely vibrational part of the Watson Hamilto-nian was expanded.

4.1. POTENTIAL ENERGY SURFACE

We used the ab initio potential energy surface ofLee, Martin, and Taylor [63] obtained by CCSD(T)singles and doubles coupled-cluster method withperturbation estimate of the effect of connected tri-ple excitations and Dunning’s correlation-consis-tent basis sets [72]. Different levels of accuracy ofthese calculations were used for the equilibriumgeometry (Dunning’s aug-cc-pVQZ basis set com-pletely uncontracted plus 1p3d2f on the carbon andcore correlation), the quadratic force constants(Dunning’s cc-pVQZ basis set), and the cubic andquartic force constants (Dunning’s cc-pVTZ basisset).

This PES is expanded in terms of internal bondand angle coordinates. We used the Mathematicasymbolic algebra package to transform analyticallythe original PES into a normal Cartesian coordinateexpansion as required for using the Watson Ham-iltonian. In doing so we have taken into account thedifferent values of the equilibrium geometry re cor-responding to the different levels of accuracy of theelectronic calculations.


254 VOL. 93, NO. 3

The coordinate transformation being a nonlinearone, we studied the effect of re-expanding the orig-inal quartic potential to different orders as Maessenand Wolfsberg did for water [73]. These authorsused a numerical re-expansion method; however,our analytic results perfectly parallel their findings.There is a large variation in ZPE of about 20 cm�1

from going to a fourth-order to a fifth-order re-expansion. Exploratory calculations using a sixth-order expansion have shown a damping of the ZPEvariation by an order of magnitude. In Table I wepresent the results of well-converged vibrational

calculations (see the following sections) corre-sponding to the fourth and fifth orders of re-expan-sion. The W50yC46 column of Table I shows thatthe wave numbers are close to those obtained witha similar VCI calculation performed with the pro-gram MULTIMODE and for the same PES but re-expanded using numerical fitting [5, 69], that is, thebending modes are close to their experimental val-ues but the stretching mode fundamentals andovertones are too large by more than 1%. As a resultthe order of the energy levels are wrongly pre-dicted, in particular for the �3 � �4 spectral lines.

TABLE I ______________________________________________________________________________________________Energies (cm�1) of the vibrational levels calculated at different orders of rotational correction andBorn–Oppenheimer potential re-expansion in terms of mass-weightened Cartesian coordinates.

Assignment Irreps. W4xyC46 W40yC46 W41yC46 W50yC46 Experimentala

�4 F2 1297.1 1308.9 1308.7 1310.9 1310.8�2 E 1519.2 1528.4 1528.2 1532.4 1533.32�4 A1 2570.5 2588.8 2588.9 2578.2 2587.32�4 F2 2585.7 2611.2 2611.0 2609.5 2614.22�4 E 2596.6 2622.3 2621.8 2628.4 2624.6�2 � �4 F2 2805.7 2827.0 2827.0 2829.5 2838.2�2 � �4 F1 2813.8 2838.6 2838.2 2844.1 2846.0�1 A1 2924.4 2925.0 2925.0 2952.9 2916.5�3 F2 3015.2 3027.0 3026.3 3059.1 3019.52�2 A1 3033.8 3051.7 3051.5 3061.5 3064.42�2 E 3036.5 3054.7 3054.4 3062.5 3065.13�4 F2 3848.0 3878.6 3878.8 3856.8 3870.53�4 A1 3872.8 3913.5 3913.3 3898.5 3909.23�4 F1 3882.3 3923.5 3923.2 3925.0 3920.53�4 F2 3894.7 3935.8 3935.4 3940.5 3930.5�2 � 2�4 E 4074.4 4104.8 4105.1 4102.5 4105.2�2 � 2�4 F1 4088.5 4127.7 4127.9 4131.2 4128.6�2 � 2�4 A1 4102.6 4139.0 4138.9 4141.3 4133.0�2 � 2�4 F2 4102.2 4142.4 4142.2 4141.5 4142.9�2 � 2�4 E 4110.6 4150.8 4150.5 4156.8 4151.2�2 � 2�4 A2 4116.5 4159.1 4158.3 4168.5 4161.9�1 � �4 F2 4213.1 4226.1 4225.8 4263.7 4223.5�3 � �4 F2 4290.8 4320.0 4318.4 4363.0 4319.2�3 � �4 E 4293.7 4326.1 4324.9 4365.1 4322.2�3 � �4 F1 4287.9 4319.9 4318.3 4364.9 4322.6�3 � �4 A1 4317.2 4328.2 4327.9 4364.5 4322.72�2 � �4 F2 4318.8 4348.2 4348.2 4353.4 4348.2�2 � �4 F1 4325.9 4360.0 4360.0 4365.9 4363.2�2 � �4 F2 4332.3 4370.0 4369.4 4378.8 4379.�1 � �2 E 4429.9 4440.5 4440.4 4473.4 4446.4�2 � �3 F1 4512.1 4539.9 4538.4 4579.9 4537.6�2 � �3 F2 4522.7 4542.5 4541.5 4584.6 4543.83�2 E 4550.9 4577.6 4577.3 4591.9 4592.03�2 A2 4556.2 4583.5 4583.0 4594.5 4595.33�2 A1 4556.1 4583.4 4582.9 4594.4 4595.6

a See reference [64].



However, inspection of the W40yC46 column ofTable I, which presents the results of the same VCIfor the fourth-order re-expanded PES, shows amuch weaker discrepancy entailing the stretchingswhereas the description of the bendings is stillgood. Therefore, in this work we stuck to limitedseries expansion rules, that is, we have used thefourth-order re-expansion in terms of Cartesiannormal coordinates of the fourth-order initial PESso as to take advantage of this error compensationeffect. This is a simpler way of dealing with thestretching mode problem than the adjusted Morserepresentation of Ref. [70] and is good enough forour present purpose of investigating the rotationalenergy levels calculated by using the generalizedperturbation theory.

4.2. EXPANSION ORDERS OF THE �MATRIX

The � matrix in the Watson Hamiltonian expres-sion is a function of the normal coordinates, whichcan be developed as a series expansion according toEq. (14). The program WATCOR implements thisformula. It requires as input the equilibrium geom-etry, the matrix relating the normal coordinates tothe mass-weighted Cartesian atomic displacements,and the coefficients of the m-order Taylor expansionof the Born–Oppenheimer potential in terms of thenormal coordinates (for a given integer m). It out-puts:

1. The coefficients of the normal coordinates and

conjugate momenta power expansion for thepurely vibrational part of the Watson Hamil-tonian (Born–Oppenheimer potential, Watsonterm, and Coriolis terms) corresponding tothe �-matrix series truncated at (an arbitrary)order n.

2. The coefficients of the normal coordinates andconjugate momenta power expansion for thevibrational part of the Watson Hamiltonian infactor of the rotational operators (�i�x, �i�y,�i�z, �x

2, �x�y � �y�x, �x�z � �z�x, �y2,

�y�z � �z�y, �z2) corresponding to the �-ma-

trix series truncated at (an arbitrary) order l.

We explained in the previous section why wehave chosen to set m � 4 for the study of thegeneralized perturbation method. Here, we estab-lish the suitable choice for the other two parame-ters. The W4nyC46 columns of Table I show theresults of the same VMFCI hierarchical scheme forthe same PES re-expansion and different choices ofn. It is clearly necessary to include the corrections ofrotational origin (essentially the Coriolis terms) atorder n � 0 because they change the wave numbersbelow 4600 cm�1 by as much as 10–40 cm�1 withrespect to the n � x calculation. However, pushingthe calculation to the next order adds 1177 moreterms to the vibrational Hamiltonian expression buthardly the energy levels affected. The most affectedwave numbers vary by less than 2 cm�1. Therefore,in the following n was set to 0.

The effect of the rotational operator expansionorder l on the rotational energy level of the vibra-

TABLE II ______________________________________________________________________________________________Energies (cm�1) of the rotational levels for increasing orders of the � matrix expansion.

ord1 ord2 ord3 ord4 ord5

J � 0 0.0 0.0 0.0 0.0 0.0J � 1 10.36256 10.48396 10.47804 10.48008 10.47993J � 2 31.08476 31.44916 31.43128 31.43742 31.43697

31.08504 31.44942 31.43155 31.43769 31.43724J � 3 62.16074 62.89015 62.85406 62.86635 62.86544

62.16185 62.89118 62.85513 62.86742 62.8665262.16324 62.89247 62.85648 62.86877 62.86786

J � 4 103.5803 104.7975 104.7365 104.7571 104.7555103.5822 104.7993 104.7384 104.7589 104.7574103.5836 104.8006 104.7398 104.7603 104.7588103.5878 104.8045 104.7438 104.7643 104.7628

The rotational Hamiltonian integrals have been calculated with the vibrational calculations labeled W40lC44, where l is the � matrixexpansion order, which include 2335 basis functions in the final contraction step. The rotational levels are calculated at the fourthorder of perturbation with summations truncated at 220 basis functions.


256 VOL. 93, NO. 3

tional ground state is presented in Table II. Only theJ � 5 energy levels are tabulated because this issufficient to determine the usual effective rotationalHamiltonian up to the octic centrifugal constants.We observe an oscillatory convergence behaviorthat is interesting to bound the truncation error.The order l � 0 is a limiting case where the rotationsare not coupled to the vibrations and where it doesnot make sense to apply the fourth-order general-ized perturbation. The spectrum corresponding tothis case would be that of the zero-order perturba-tion spectrum presented later in Table IV. For sym-metry reasons, the corrections obtained at order l �1 are exclusively due to the totally symmetrical A1vibrational mode. Consequently, they can only pro-duce an effective rotational Hamiltonian that is apolynomial in �2. At order l � 2 this symmetryrestriction is lifted and we see that 4 digits areconverged for the J � 1 level and three digits for theother levels. Going to higher orders, the fluctua-tions have less amplitude and it seems that 5 digitsare consistently converged at order l � 4 exceptperhaps for the highest J � 3 level (but, order 6should give a larger value than that of order 5).Because the ab initio value of the equilibrium bondlength is only given to 4 digits, we can reasonablylimit the calculations to order l � 4 in the following.

Now, we have fixed the form of the WatsonHamiltonian expansion that we will use. In the nextsection, we explain how we constructed the zero-order vibrational wave functions.

4.3. CONTRACTION SCHEME FOR THEVIBRATIONAL CONFIGURATIONINTERACTION

The VMFCI procedure described in Section 2 hasbeen applied to the Watson Hamiltonian expan-sions detailed in the previous sections. The moreefficient contraction scheme has been determinedfrom test calculations, to be detailed elsewhere [50].In the present work we used the more naturalscheme illustrated in Figure 1. It consists of sixsuccessive VMFCI steps. The components of thedegenerate modes are first contracted to preparesymmetry-adapted eigenvectors to be used in thenext steps. The initial set of 9 modals is thus re-duced to 4 sets of product functions describing thenormal modes of methane, respectively, of A1, E, F2,and F2 symmetries (step 0). Then, this first VMFCIis repeated 3 times (steps 1 to 3) to perform a“normal-mode SCF calculation.” This achieves aconvergence for the fundamentals to about 1 cm�1.

Step 4 contracts independently the stretchingmodes (�1 and �3) and the bending modes (�2 and�4), while the last step contracts stretching andbending modes in the final n-mode VCI calculation.

In line with Section 2.3, we can write the abovecontraction scheme in the compact form: {{{1}4, {4, 5,6}4}, {{2, 3}4, {7, 8, 9}4}}, where {1} corresponds to theA1 normal mode, {2, 3} to the E normal mode, and{4, 5, 6} and {7, 8, 9} to the F2 normal modes. In thisnotation, the nested brackets represent the succes-sive contractions, the first steps being described inthe innermost brackets. The exponents indicate thenumber of times a given VMFCI step is iterated.

Only two types of vibrational calculations havebeen used in this work to investigate the predictionof vibrational and rotational spectroscopic proper-ties by our method. They share the contractionscheme of Figure 1 and the number of harmonicoscillator functions in the initial basis set for thestretching modals, that is, 7. They differ by thenumber of harmonic oscillators in the initial basisset for the bending modals and the selection of theproduct functions to be used in the two last con-traction steps.

The first type corresponds to the calculationslabeled WmnlC46. These calculations start with 9harmonic oscillator functions per bending modal.In step 4, the stretching product functions retainedare those made of normal-mode components whosesum of the �1 and �3 mode wave numbers is lessthan 24,000 cm�1. The bending product functionsretained are those made of normal-mode compo-nents that satisfy the following two criteria: (1) Thesum of the �2 and �4 modes wave numbers is lessthan 14,000 cm�1 and (2) the �4 mode wave numberis less than 12,000 cm�1 with respect to the ZPE ofthe mode. In step 5, the product functions selectedare those made of VMFCI components satisfying (1)

FIGURE 1. Contraction scheme used for this work.



the sum of stretching VMFCI wave number and thebending VMFCI wave number is less than 12,300.0cm�1 and (2) the bending VMFCI wave number is lessthan 9300 cm�1 with respect to its ZPE. All thresholdswere chosen so as to retain in the final calculation allthe states corresponding to an approximate stretchingquantum number of 4 and an approximate bendingquantum number of 6 while keeping the VCI matrixto diagonalize as small as possible. The size of thefinal VCI matrix was 4025 � 4025, and its 1196 lowesteigenvalues (counted with their degeneracy) corre-sponded to approximate quantum numbers withinour criterion.

The second type corresponds to the calculationslabeled WmnlC44. These calculations start with 7 har-monic oscillator functions per bending modal andhave more stringent truncation thresholds becausethey were designed to retain only states with an ap-proximate bending quantum number of 4. The thresh-olds were, respectively, 22,000, 9200, and 8050 cm�1

for the sum of stretching normal modes wave num-bers, the sum of bending normal modes wave num-bers, and the �4 mode wave number in step 4, respec-tively, and 12,200 and 6200 cm�1 for, respectively, thesum of stretching VMFCI plus bending VMFCI wavenumber and the bending VMFCI wave number. Thesize of the final VCI matrix was 2335 � 2335.

The vibrational calculations were performedwith the code CONVIV. This program takes asinput the coefficients of the normal coordinates andmomenta operators, Eq. (12), punched by the pro-gram WATCOR and the contraction/truncation

scheme designed by the user. It outputs the vibra-tional eigenvalues and eigenfunctions at eachVMFCI step and punches an integral file needed byROTEFF, the next program in our chain. ROTEFFimplements the perturbation formulas of Section3.2 for a given perturbation order (less than 4) anda given vibrational state chosen by the user. It out-puts the rotational effective Hamiltonian, which isin turn diagonalized in a spherical harmonic basisset using a Mathematica script called MakeSpec.This last step, which gives the rotational energylevels, is cheap in terms of computational require-ments because the matrices of the rotational angu-lar momentum operators have been calculated oncefor all and stored up to J � 18.

Remark 4. Our program is not designed to calcu-late rotational constants that are not always unambig-uously defined [74] but energy levels and transitionmoments (when a dipole moment surface is avail-able). This is why we have not presented tables ofrotational and centrifugal distortion constants. How-ever, these constants can be derived from our energylevel tables by spectroscopy fitting packages.

4.4. TRUNCATION THRESHOLD FOR THESUMS OVER VIBRATIONAL FUNCTIONS

Table III shows the convergence of the rotationalenergy levels of the vibrational ground state withincreasing truncation thresholds for the summa-tions over vibrational states in the perturbation se-

TABLE III _____________________________________________________________________________________________Energies (cm�1) of the rotational levels calculated at the fourth order of perturbation with different seriestruncation thresholds.

0 cm�1 (1f) 3800 cm�1 (25f) 6400 cm�1 (220f) 7700 cm�1 (536f) 9200 cm�1 (1196f)

J � 0 0.0 0.0 0.0 0.0 0.0J � 1 10.49362 10.49311 10.48009 10.47900 10.47899J � 2 31.48087 31.47651 31.43743 31.43419 31.43415

31.48087 31.47677 31.43770 31.43446 31.43442J � 3 62.96174 62.94455 62.86637 62.85992 62.85985

62.96174 62.94563 62.86744 62.86100 62.8609362.96174 62.94697 62.86878 62.86235 62.86228

J � 4 104.9362 104.8875 104.7571 104.7464 104.7463104.9362 104.8893 104.7590 104.7483 104.7482104.9362 104.8907 104.7603 104.7497 104.7495104.9362 104.8947 104.7643 104.7537 104.7536

The energy thresholds are expressed in cm�1 and the numbers of eigenfunctions corresponding to eigenvalues less or equal to thethresholds are given in parentheses. The rotational Hamiltonian integrals have been calculated with the vibrational calculationreferred to in the text as W404C46, which includes 4025 basis functions in the final contraction step.


258 VOL. 93, NO. 3

ries. The vibrational functions used to build theeffective rotational Hamiltonian were obtained bydiagonalizing in the final contraction step a 4025 �4025 VCI matrix that includes all vibrational statesbelow 9200 cm�1. The truncation thresholds of thetable were chosen to fall in the gaps of the vibra-tional spectrum, that is, the columns of Table IIIcorrespond to summations involving all stretchingand bending modes with approximate quantumnumber nbend � 1 (limiting case corresponding tofirst-order perturbation because the vibrationalground state is excluded from the summation inhigher-order perturbative contributions), nbend � 3,nbend � 5, nbend � 6, nbend � 7, respectively. Notethat an nbend � 2 � k � 1 truncation corresponds toan nstretch � k � 1 truncation in this energy range.

We observe a general lowering of the levels asthe number of functions used in the summationincreases. This lowering grows with J approxi-mately as J( J � 1). Therefore, this is to be related toa decrease of the rotational constant B in factor ofJ( J � 1) in the usual effective rotational Hamilto-nian. The overall convergence is satisfactory be-cause 6 digits seem converged with a summationrestricted to 535 functions (the ground state beingexcluded) except for 2 sublevels of the J � 3 level,where only 5 digits are converged. The convergenceof the sublevel splittings is even more encouragingbecause they are almost converged with only 25functions. This means that the centrifugal distortionconstants of the corresponding effective Hamilto-nian are rapidly converged. Consequently, a moredrastic truncation can be considered for the sum-

mations in the third- and fourth-order perturbativecontributions because, as we shall see in the nextsection, they have little influence on the rotationalconstant. This is a particularly important featurebecause the computational effort for the nth-orderperturbation grows like the number of functionsinvolved in the summation to the power n � 1.

4.5. CONVERGENCE WITH INCREASINGPERTURBATION ORDERS

The convergence of the rotational energy levelsof the vibrational ground state with increasing per-turbation order is displayed in Table IV. The zeroorder corresponds here to the spherical rotor ob-tained for the equilibrium rotational constant Be

and the first order to the spherical rotor obtained byaveraging the � matrix diagonal element (times afactor 1/2) over the vibrational ground state. It isclear that the averaging has a large effect and im-proves drastically the rotational levels. However,the degeneracy is not lifted at this level becausequartic centrifugal distortion only appears in thesecond-order corrective term Eq. (36), as alreadynoted in Section 3.2. As a matter of fact, it is neces-sary that at least two vibrationally averaged H1 bemultiplied to produce a term of fourth order in theangular momentum components, ��, in our effec-tive rotational Hamiltonian. It also proves sufficienthere as the degeneracy is removed at order 2. Be-sides, we observe that the levels are converged to 5digits (almost 6) at the second order. The thirdorder, which introduces sextic centrifugal distor-

TABLE IV _____________________________________________________________________________________________Energies (cm�1) of the rotational levels for increasing perturbation orders.

ord0 ord1 ord2 ord3 ord4

J � 0 0.0 0.0 0.0 0.0 0.0J � 1 10.63296 10.49362 10.48010 10.48008 10.48008J � 2 31.89887 31.48087 31.43746 31.43742 31.43742

31.89887 31.48087 31.43772 31.43769 31.43769J � 3 63.79775 62.96174 62.86645 62.86637 62.86635

63.79775 62.96174 62.86749 62.86745 62.8674263.79775 62.96174 62.86879 62.86879 62.86877

J � 4 106.3296 104.9362 104.7573 104.7571 104.7571106.3296 104.9362 104.7591 104.7590 104.7589106.3296 104.9362 104.7604 104.7603 104.7603106.3296 104.9362 104.7643 104.7644 104.7643

The rotational Hamiltonian integrals have been calculated with the vibrational calculations labeled W404C44, which include 2335basis functions in the final contraction step. However, the summations in the perturbative expansion were truncated at 220 basisfunctions. ord0 corresponds to the spherical rotor obtained for the equilibrium rotational constant Be.



TABLE V ______________________________________________________________________________________________Energies (cm�1) with respect to the J � 0 level given by the STDS package, which consistently reproducesthe transitions tabulated in the HITRAN database [77] with 7 or more accurate digits (first column) andpredicted with our method (second column).

J Irreps. STDS This work J Irreps. STDS This work

1 F1 10.481648 10.481652 E 31.442121 31.442132 F2 31.442387 31.442403 F1 62.875779 62.875813 F2 62.876841 62.876893 A2 62.878169 62.878244 A1 104.772844 104.77294 F1 104.774704 104.77484 E 104.776033 104.77614 F2 104.780011 104.78025 F1 157.124338 157.12445 F2 157.127927 157.12815 E 157.137194 157.13755 F1 157.138918 157.13936 E 219.913464 219.91366 F2 219.915049 219.91526 A2 219.919854 219.92016 F2 219.936765 219.93736 F1 219.941261 219.94196 A1 219.945233 219.94597 F1 293.122987 293.12327 F2 293.126548 293.12687 A2 293.154197 293.15497 F2 293.164567 293.16557 E 293.170133 293.17117 F1 293.178678 293.17988 A1 376.730436 376.73078 F1 376.733719 376.73418 E 376.735645 376.73608 F2 376.785873 376.78708 F1 376.804777 376.80628 E 376.821288 376.82308 F2 376.826268 376.82809 F1 470.716958 470.71739 F2 470.720342 470.72089 E 470.798974 470.80079 F1 470.805281 470.80719 A1 470.830962 470.83319 F1 470.855000 470.85759 F2 470.865058 470.86779 A2 470.872917 470.8757

10 E 575.051266 575.051710 F2 575.052656 575.053110 A2 575.055668 575.056110 F2 575.170076 575.172410 F1 575.184305 575.186810 A1 575.222922 575.2261

10 F1 575.259785 575.263410 E 575.271915 575.275810 F2 575.285416 575.289511 F1 689.705386 689.705811 F2 689.707734 689.708211 A2 689.862325 689.865211 F2 689.876914 689.880011 E 689.886292 689.889611 F1 689.956854 689.961311 F2 690.017637 690.022911 E 690.039743 690.045411 F1 690.049396 690.055212 A1 814.646210 814.646512 F1 814.648073 814.648412 E 814.649035 814.649412 F2 814.866671 814.870512 F1 814.884499 814.888612 E 814.993340 814.999212 F2 815.008253 815.014312 A2 815.089170 815.096412 F2 815.116077 815.123712 F1 815.131686 815.139612 A1 815.143782 815.151913 F1 949.841648 949.841813 F2 949.843144 949.843313 E 950.129491 950.134213 F1 950.136556 950.141413 A1 950.153920 950.159013 F1 950.304979 950.312513 F2 950.337020 950.345013 A2 950.385253 950.394013 F2 950.487439 950.497613 E 950.505128 950.515613 F1 950.522428 950.533214 E 1095.251425 1095.251214 F2 1095.251995 1095.251714 A2 1095.253149 1095.252914 F2 1095.619349 1095.625014 F1 1095.632143 1095.638014 A1 1095.828810 1095.837914 F1 1095.869302 1095.878914 E 1095.894821 1095.904814 F2 1095.981538 1095.993014 F1 1096.133087 1096.146614 E 1096.157058 1096.171014 F2 1096.171074 1096.185215 F1 1250.836389 1250.8356

(continued)


260 VOL. 93, NO. 3

tion, ensures the convergence of the sixth digit withrespect to the fourth order, which contains octiccentrifugal distortion and serves as a reference.

As already mentioned the ab initio equilibriumgeometry is at best 4 digits accurate. Therefore, fora purely ab initio prediction one could stop at thesecond order of perturbation. However, as we shallsee in the next section, it is possible to achieve abetter accuracy by means of a simple one-parame-ter scaling that implies the use of at least third-order perturbation for error consistency reasons.

4.6. COMPARISON WITH EXPERIMENT

4.6.1. Vibrational Energy Levels

All the energy levels up to 9200 cm�1 have beencomputed at different levels of accuracy and are

available upon request. Those below 4600 cm�1 aretabulated in Table I.

No systematic comparison with previous the-oretical work on the vibrational energy levels ofmethane will be given in the present work, whichis essentially devoted to the rotational part. Sucha comparison merits a throughout analysis of ourvibrational contraction approach and a carefulexamination of the PESs used by the differentauthors. A detailed analysis of the vibrationalresults will be proposed in a forthcoming article[50].

As explained in Section 4.1 our vibrational en-ergy levels were not meant to have state-of-the-artaccuracy. However, they come out to be fairly sat-isfactory. The calculations with the PES re-ex-panded at fourth order in rectilinear normal coor-

TABLE V ______________________________________________________________________________________________(Continued)

J Irreps. STDS This work J Irreps. STDS This work

15 F1 1250.836389 1250.835615 F2 1250.837251 1250.836415 A2 1251.291502 1251.298015 F2 1251.302087 1251.308715 E 1251.307878 1251.314615 F1 1251.590702 1251.602115 F2 1251.646861 1251.659015 E 1251.779934 1251.794215 F1 1251.807331 1251.822015 A1 1252.002106 1252.019515 F1 1252.027835 1252.045715 F2 1252.046484 1252.064715 A2 1252.061680 1252.080116 A1 1416.553666 1416.552016 F1 1416.554305 1416.552616 E 1416.554626 1416.552916 F2 1417.119957 1417.127416 F1 1417.129360 1417.136916 E 1417.499583 1417.513216 F2 1417.519829 1417.533716 A2 1417.579350 1417.594016 F2 1417.753110 1417.770716 F1 1417.807375 1417.825816 A1 1417.865247 1417.884616 F1 1418.099134 1418.121716 E 1418.118763 1418.141716 F2 1418.137412 1418.160617 F1 1592.359243 1592.356317 F2 1592.359714 1592.3568

STDS [75] makes use of the rotational constants and the centrifugal distortion constants (up to the octic) determined by fittingexperimental spectra. Our estimation is purely ab initio except for the scaling of �0��0� (which can be seen roughly as scaling theequilibrium length by an empirical factor of 0.99987).

17 E 1593.044952 1593.053117 F1 1593.048625 1593.056817 A1 1593.056319 1593.064617 F1 1593.523300 1593.539217 F2 1593.562048 1593.578417 A2 1593.784741 1593.804917 F2 1593.878631 1593.900117 E 1593.928421 1593.950617 F1 1594.026232 1594.050117 F2 1594.344986 1594.373317 E 1594.367216 1594.395917 F1 1594.383307 1594.412318 E 1778.204858 1778.200318 F2 1778.205029 1778.200418 A2 1778.205371 1778.200818 F2 1779.024766 1779.033318 F1 1779.030621 1779.039318 A1 1779.597416 1779.615318 F1 1779.631419 1779.649818 E 1779.651450 1779.670018 F2 1779.984723 1780.008918 F1 1780.122654 1780.148718 E 1780.269506 1780.298018 F2 1780.312383 1780.341618 A2 1780.707364 1780.742118 F2 1780.729385 1780.764618 F1 1780.747755 1780.783318 A1 1780.763728 1780.7995



dinates predict correctly the order of the energylevels with the following exceptions:

▪ The �2 � 2�4 A1 and F2 are swapped when norotational corrections are taken into account.

▪ The �3 � �4 F1 state is consistently foundbelow the �3 � �4 F2 and not above the �3 � �4E state as in the “Experimental” column ofTable I.

▪ The 3�2 A1 state is systematically found belowthe 3�2 A2 compared to the “Experimental”column of Table I.

The calculation with the PES re-expanded at fifthorder in rectilinear normal coordinates is less reliablebecause of the stretching mode discrepancy; however,the �2 fundamental and overtones are in better agree-ment with the experiment than the other calculationsand it, too, predicts the 3�2 A1 state below the 3�2 A2.

Wang and Sibert [66], who used the same PES aswe did but re-expanded in a different way to im-prove the stretching mode description, have notreported the inversions listed above. Schwenke andPartridge [71] did observe the 3�2 inversion withtheir quartic force field, but the order derived fromexperiment is retrieved with their octic force field.The small energy differences involved are well be-low our error bar and prevent any definite conclu-sion from our calculations, but in view of the con-sistency of the other overtone predictions we wouldlike to draw attention to possible misassignments inthe 3�2 and/or �3 � �4 bands.

4.6.2. Rotational Energy Levels of theVibrational Fundamental

Our predicted rotational energy levels up to J � 18are tabulated in Table V together with those predictedby the STDS package of Champion and Wenger [75].This code uses the effective rotational and centrifugalconstants refined on experimental spectra [76] to gen-erate the rotational energy levels. It reproduce consis-tently the observed transitions tabulated in the HIT-RAN database [77] with 7 or more accurate digits andthe energy levels themselves are expected to have an8-digit accuracy. Because we cannot expect to achievesuch an accuracy with our calculations we have con-veniently taken the STDS numbers for the “experi-mental” references.

With the equilibrium distance value re �1.0862 � 0.0005 Å and the force field of Ref. [63],our perturbation method proved able to reproduce

the STDS numbers with 4 accurate digits. Therefore,the step limiting the accuracy of a purely ab initioprediction is not the rovibrational treatment but theelectronic calculation (and possibly the Born–Op-penheimer approximation). We tried to circumventthis limitation by scaling the expectation value ofthe �-matrix diagonal element �0��0� (with a scal-ing factor of roughly 1.0002535485) so as to obtainthe J � 1 level with 7 accurate digits. It is the levelsscaled in this way that are reported in Table V.Their unscaled counterparts were calculated at thefourth order of perturbation with � expanded at thefourth order and summations truncated at 1195excited vibrational eigenfunctions. Our attemptproves successful because 6 accurate digits are ob-tained in the lower part of the spectrum and 5 digitsare retained up to the highest part with relative errorsless than 2 � 10�5. Thus, it is possible to gain twodecimal places with respect to the pure ab initio pre-cision by simply adjusting one parameter. By makingthe assumption that scaling �0��0� amounts in firstapproximation to scaling its constant part ��

e , whichin turn amounts to scaling re by the square root of theinverse of the scaling factor, we derived an equilib-rium distance re � 1.08606 Å within the error bar ofthe ab initio predicted value.

5. Conclusion

We implemented a new method for the calcula-tion of molecular rotation–vibration spectra. Themethod is immediately applicable to a wide rangeof molecules, whether they are spherical, symmet-rical, or asymmetrical top systems. It is not limitedto small values of the rotational quantum numbers,the computational effort for the rotational levelsbeing almost independent of J. The method hasbeen applied to methane and the results have con-firmed its great potential interest.

The VMFCI technics permits useful and flexibletruncations of the vibrational basis set that reducethe size of the vibrational problem. In methane ithas been demonstrated [50] that all the intermodecouplings can be accounted for in a satisfactorymanner by contracting at most two normal modesat a time. For example, the five-mode couplinginvoked by Carter and Bowman [69] to explainsome of their results does not prevent the conver-gence of the affected vibrational energy levels inour calculation. The origin of this success is theinclusion in the VMFCI equations of the mean-fieldeffect of the coupling with the other vibrational


262 VOL. 93, NO. 3

degrees of freedom, which prepares compact basissets of good quality for the final contraction step. Inthe light of this observation it seems reasonable tohope that larger molecules than a penta-atomic sys-tem will be amenable to a VMFCI treatment thatconverges the vibrational energies to the spectro-scopic accuracy, so that it will be the PES qualityalone (not the vibrational treatment) that will limitthe precision of the ab initio predictions.

The study of the convergence of our perturbationscheme have shown that:

1. The �-matrix need to be expanded to fourthorder to converge five digits. This is a largerorder than what is usually considered in thetraditional perturbation theory. At this order,the expansion of H1 contains 10,838 terms infactor of the 9 rotational operators of Section 4.2.

2. The number of vibrational eigenfunctions tosum over in the perturbation series can belimited to 535 in the methane case. This cor-responds to a difference in approximate bend-ing quantum number of 5 with respect to theground state. However, this number can bereduced in the third- and higher-order pertur-bative corrections.

3. The second order of perturbation is sufficientwith respect to the accuracy of the electroniccalculations, the third order can be useful inconjunction with an empirical scaling of theequilibrium geometry, and the fourth order ofperturbation will not be justified until theforce field accuracy has improved.

4. The predictions of the rotational energy levelsare good even for energy levels well above thelowest excited vibrational states, that is, whenthe initial hypothesis justifying our perturba-tive expansion is not met.

There are two computational bottlenecks atpresent in our approach:

1. A memory space bottleneck when the diagonal-ization of a large VCI matrix is needed. As amatter of fact, we need a large number of vibra-tional eigenfunctions to sum over in the pertur-bation series. Many of them can be degenerateor quasidegenerate. We have used so far thetred2 and tql2 routines of the Lapack fortranlibrary, which calculate and provide all the eig-envalues and eigenvectors of a given matrix.Less memory space-consuming diagonalizers

based on the Davidson algorithm [78] have beentried without success in our problem.

2. A CPU time bottleneck in the last contractionstep when computing and summing the inte-grals over vibrational eigenfunctions requiredin the perturbation energy expression. In theW404C46 run we consumed an average of0.03858 s of CPU time per pair of eigenfunctionson a Compaq digital EV6.7 alpha CPU (700MHz, 8 Mb cache memory). That is 7 h and 40min of CPU time for 1196 functions and 1.5 h for536 functions. There were 10,876 different typesof integrals to be calculated and summed for the454 terms of the H0 expansion and the 10,838terms of the H1 expansion.

The extension of the present study to the otherisotopomers of methane and excited vibrationalstates is in progress. It will allow us to assess thetransferability of the empirical scale factor obtainedhere and the generality of the conclusions reachedin this study. The computation of transition mo-ments is also planned. We hope that this will lead toa better understanding of the methane spectroscopyin the near future and of other, possibly larger,systems in the long run.

ACKNOWLEDGMENTS

The authors thank the CGRI-CNRS-FNRS for a“Tournesol” French–Belgium cooperation grantand also G. Destree of the ULB/VUB computercenter for his help in some computational prob-lems. J.L. acknowledges the Fonds National de laRecherche Scientifique for financial support (FRFCand IISN contracts).

References

1. Harding, L. B.; Ermler, W. C. J Comput Chem 1985, 6, 13.2. Gaw, J. F.; Willets, A.; Green, W. H.; Handy, N. C. In:

Bowman, J., ed. Advances in Molecular Vibrations and Col-lision Dynamics; JAI Press; Greenwich, CT, 1991.

3. Munoz-Caro, C.; Nino, A. QCPE Program 665.4. Dunn, K. M.; Boggs, J. E.; Pulay, P. J Chem Phys 1986, 85, 5838.5. Carter, S.; Bowman, J. M. In: Law, M. M.; Atkinson, I. A.;

Hutson, J. M., eds. Rovibrational Bound States in PolyatomicMolecules; (CCP6: Daresbury, 1999); pp 60–67.

6. Demaison, J.; Nemes, L. J Mol Spectrosc 1979, 55, 295.7. Hollenstein, H.; Marquardt, R. R.; Quack, M.; Suhm, M. A.

J Chem Phys 1994, 101, 3588.8. Watson, J. K. G. Mol Phys 1970, 19, 465–487.



9. Watson, J. K. G. Mol Phys 1968, 15, 479–490.10. Sugny, D.; Joyeux, M. J Chem Phys 2000, 112, 31.11. Lauvergnat, D.; Nauts, A.; Justum, Y.; Chapuisat, X. J Chem

Phys 2001, 114, 6592.12. Searles, D.; Von Nagy-Felsobuki, E. Ab Initio Variational

Calculations of Molecular Vibrational–Rotational Spectra,Lecture Notes in Chemistry; Springer-Verlag: Berlin, 1993.

13. Herman, M.; Lievin, J.; Vander Auwera, J.; Campargue, A.In: Prigogine, I.; Rice, S. A., eds. Advances in ChemicalPhysics 108; Wiley: New York, 1999, pp 1–431.

14. Gatti, F.; Lung, C.; Leforestier, C.; Chapuisat, X. J Chem Phys1999, 111, 7236.

15. Viel, A.; Leforestier, C. J Chem Phys 2000, 112, 1212.16. Van Vleck, J. H. Phys Rev 1929, 33, 467.17. Nielsen, H. H. Phys Rev 1941, 60, 794.18. Mills, I. M. In: Narahari Rao, K.; Mathews, C. W., eds.

Molecular Spectroscopy: Modern Research; Academic Press:New York, 1972; pp 115–140.

19. Aliev, M. R.; Watson, J. K. G. In: Narahari Rao, K., ed.Molecular Spectroscopy: Modern Research, Vol. 3; AcademicPress: Orlando, FL, 1985, pp 1–67.

20. Shaffer, W. H.; Nielsen, H. H. Phys Rev 1939, 56, 188.21. Carney, G. D.; Sprandel, L. L.; Kern, C. W. In: Prigogine, I.;

Rice, S. A., eds. Advances in Chemical Physics, Vol. 37;Wiley: New York, 1978; pp 305–379.

22. Bryukhanov, V. N.; Makushkin, Y. S. Opt Spectrosc 1974, 36,469–474.

23. Makushkin, Y. S. Opt Spectrosc 1974, 37, 662–667.24. Voitsekhovskaya, O. K.; Makushkin, Y. S. Opt Spectrosc

1975, 39, 32–37.25. Maessen, B.; Wolfsberg, M. J Phys 1985, 89, 3876.26. Brodersen, S.; Lolck, J.-E. J Mol Spectrosc 1987, 126, 405–426.27. Pouchan, C.; Zaki, K. J Chem Phys 1997, 107, 342.28. Born, M.; Oppenheimer, J. R. Ann Phys 1927, 84, 457.29. Pickett, H. M. J Chem Phys 1972, 56, 1715.30. Calogeracos, A.; Castillejo, L. Mol Phys 1993, 80, 1359.31. Littlejohn, R. G.; Reinsch, M. Rev Mod Phys 1997, 69, 213.32. Louck, J. D. J Mol Spectrosc 1976, 61, 107.33. Darling, B. T.; Dennison, D. M. Phys Rev 1940, 57, 128.34. Podolsky, B. Phys Rev 1928, 32, 812.35. Wilson, E. B. Jr.; Howard, J. B. J Chem Phys 1936, 4, 260.36. Whitehead, R. J.; Handy, N. C. J Mol Spectrosc 1975, 55, 356.37. Romanowski, H.; Bowman, J. M.; Harding, L. B. J Chem Phys

1985, 82, 4155.38. Cassam-Chenaı, P.; Pauzat, F.; Ellinger, Y. J Mol Struct Theo-

chem 1995, 330, 167.39. Gerber, R. B.; Ratner, M. A. Chem Phys Lett 1979, 68, 195.40. Roothaan, C. C. J. Rev Mod Phys 1951, 23, 69.41. Thompson, T. C.; Truhlar, D. G. J Chem Phys 1982, 77, 3031.42. Zuniga, J.; Bastida, A.; Requena, A. J Phys Chem 1992, 96, 9691.43. Schwenke, D. W. J Chem Phys 1992, 96, 3426.44. Culot, F.; Lievin, J. Theor Chim Acta 1994, 89, 227.45. Culot, F.; Laruelle, F.; Lievin, J. Theor Chim Acta 1995, 92, 211.46. Ruedenberg, K.; Schmidt, M. W.; Gilbert, M. M.; Elbert, S. T.

Chem Phys 1982, 71, 41.

47. Roos, B. O. Adv Chem Phys 1987, LXVII, 399.

48. Chen, C. L.; Maessen, B.; Wolfsberg, M. J Chem Phys 1985,83, 1795.

49. Carter, S.; Handy, N. C. Comp Phys Rep 1986, 5, 15.

50. Cassam-Chenaı, P.; Lievin, J. To be published.

51. Clabo, D. A. Jr.; Allen, W. D.; Remington, R. B.; Yamaguchi,Y.; Schaeffer, H. F. III. Chem Phys 1988, 123, 187.

52. Primas, H. Rev Mod Phys 1963, 35, 710.

53. Sibert, E. L. III. J Chem Phys 1988, 88, 4378.

54. Lukka, T. J.; Kauppi, E. J Chem Phys 1995, 103, 6586.

55. Oka, T. J Chem Phys 1967, 47, 5410.

56. Messiah, A. Mecanique Quantique Tome 2; Dunod: Paris,1964; Section XVI–17.

57. Eyring, H.; Walter, J.; Kimball, D. E. Quantum Mechanics;Wiley: New York, 1947; Chapter 14, Section 2.

58. Norris, L. S.; Ratner, M. A.; Roitberg, A. E.; Gerber, R. B.J Chem Phys 1996, 105, 11261.

59. Survenev, A. A.; Goodson, D. Z. Chem Phys Lett 1997, 269,177, and references therein.

60. Dinelli, B. M.; Miller, S.; Achilleos, N.; Lam, H. A.; Cahill, M.;Tennyson, J.; Jagod, M.-F.; Oka, T.; Hilico, J.-C.; Geballe, T. R.Icarus 1997, 126, 107.

61. Allard, F.; Hauschildt, P. H.; Baraffe, I.; Chabrier, G. Astro-phys J 1996, 465, L123.

62. Geballe, T. R.; Kulkarni, S. R.; Woodward, C. E.; Sloan, G. C.Astrophys J 1996, 467, L101.

63. Lee, T. J.; Martin, J. M. L.; Taylor, P. R. J Chem Phys 1995,102, 254.

64. Ref. [48] to Hilico’s private communication. In: Halonen, L.J Chem Phys 1997, 106, 831.

65. Wiesenfeld, L. J Mol Spectrosc 1997, 184, 277.

66. Wang, X.-G.; Sibert, E. L. III. J Chem Phys 1999, 111, 4510.

67. Venuti, E.; Halonen, L.; Della Valle, R. G. J Chem Phys 1999,110, 7339.

68. Tennyson, J.; Xie, J. In: Law, M. M.; Atkinson, I. A.; Hutson,J. M., eds. Rovibrational Bound States in Polyatomic Mole-cules; CCP6: Daresbury, 1999; pp 5–12.

69. Carter, S.; Shnider, H. M.; Bowman, J. M. J Chem Phys 1999,110, 8417.

70. Carter, S.; Bowman, J. M. J Phys Chem 2000, 104, 2355.

71. Schwenke, D. W.; Partridge, H. Spectrosc Acta A 2001, 57, 887.

72. Dunning, T. H. J Chem Phys 1989, 90, 1007.

73. Maessen, B.; Wolfsberg, M. J Phys Chem 1985, 89, 3324.

74. Perevalov, V. I.; Tyuterev, V. G.; Zhilinskii, B. I. Chem PhysLett 1984, 104, 455.

75. Champion, J. P.; Wenger, C. Spherical top data system, July 1998version, freeware available from http://www.u-bourgogne.fr/LPUB/sTDS.html; Laboratoire de Physique de l’Universitede Bourgogne Unite associee au C.N.R.S.: Cedex, France.

76. Roche, C.; Champion, J. P. Can J Phys 1991, 69, 40.

77. Rothman, L. S.; Schroeder, J.; McCann, A.; Gamache, R. R.;Wattson, R. B.; Flaud, J.-M.; Perrin, A.; Dana, V.; Mandin,J.-Y.; Goldman, A.; Massie, S.; Varanasi, P.; Yoshino, K. JQuant Spectrosc Radiat Transfer 1998, 60, 665.

78. Davidson, E. R. J Comput Chem 1975, 17, 87.


264 VOL. 93, NO. 3

Documents

Alternative perturbation method for the molecular vibration–rotation problem