A numerical approach to the internal large amplitude motion Hamiltonian of a polyatomic molecule

JOURNAL

OF MOLECULAR

SPECTROSCOPY

143, l37- 159 (1990)

A Numerical Approach to the Internal Large Amplitude Motion Hamiltonian of a Polyatomic Molecule J. PYKA, I. FOLTYNOWICZ,

AND J. MAKAREWICZ

Faculty of Chemistry. A. Mickiewicz University. PL 60-780 Poznan, Poland A fully numerical approach to the vibrational Hamiltonian of a polyatomic molecule is presented and applied to four-membered ring molecules. Large amplitude motions (LAMS) and small amplitude vibrations (SAVs) are described by curvilinear coordinates. On the grounds of the adiabatic separation of LAMS and SAVs, a very efficient method of solving the inverse spectral problem is proposed. Coordinate transformations which eliminate the kinetic LAM-SAV interaction and ensure high accuracy of the adiabatic method are discussed. o 1990 Academic press. IX. I. INTRODUCTION

Large amplitude motions (LAMS) in molecules are usually of low frequency and they are approximately separated from the high-frequency small amplitude vibrations (SAVs). Such a separability justifies the application of approximate Hamiltonians describing only LAMS. For years, simple model Hamiltonians have been applied to various kinds of LAMS such as internal rotation (I ), bending (2-6), inversional ( 710), and ring-puckering motions; for numerous examples of ring-puckering models the reader is referred to reviews (1 I, 12) and references therein. Recently, it was demonstrated that one-dimensional ring-puckering model Hamiltonians were inadequate for describing the far-infrared spectra of some ring molecules ( 13-15). Consequently, two- and three-dimensional models, taking into account selected vibrational modes treated as LAMS, have been developed by Laane and coworkers (13-18). A general theory describing the interaction of one LAM with all SAVs and overall rotation has been proposed for triatomic molecules by Hougen, Bunker, and Johns (HBJ ) (3)) and further developed by Bunker and co-workers ( 19-24) and Jensen (25). This theory has been extended to the ammonia molecule by Spirko et al. (2628). In HBJ theory, rectilinear Cartesian coordinates are used for the SAVs and a curvilinear coordinate is used for the LAM. The LAM-SAV interaction in the kinetic energy is eliminated, in a zeroth-order approximation, by using the Sayvetz (29) condition. A concurrent theory of LAM-SAV interaction has been developed by Quade (30) and Cress and Quade (31). In this theory, curvilinear coordinates are used not only for the LAM but also for the SAVs. To eliminate the LAM-SAV kinetic interaction, these coordinates are transformed directly to the new ones. The HBJ and Quade theories are difficult to apply to large polyatomic molecules of low symmetry, since it is very difficult to perform the coordinate transformation 137

0022-2852190 $3.00 Copyri~t

0

1990 by Academc

h,

Inc.

All rightsof reproduction in any form mewed

138

PYKA,

FOLTYNOWICZ,

AND

MAKAREWICZ

analytically. Also, it is practically impossible to calculate analytically the kinetic tensor G for a complicated LAM mode. In this paper we present an alternative, numerical approach to the pure internal Hamiltonian (the overall rotation is not considered). This approach can be applied to an arbitrary polyatomic molecule with several LAMS. However, curvilinear internal coordinates must be defined individually for each molecule, in a separate subroutine. The effective LAM Hamiltonian that takes into account interactions with SAVs is derived by applying the adiabatic theory (Section II). On the grounds of the adiabatic separation of LAMS from SAVs, a very efficient method of solving the inverse spectral problem is proposed in Section III. The coordinate transformations, which eliminate the kinetic LAM-SAV interaction and ensure high accuracy of the adiabatic method, are discussed in Section IV. In Section V the presented theory is applied to the disilacyclobutane molecule. In the Appendix, a definition of curvilinear coordinates for four-membered ring molecules is given. II. ADIABATIC

SEPARATION

OF LAMS FROM

SAVS

The effective Hamiltonian for LAMS can be derived by applying the adiabatic method (32-43). The adiabatic wavefunction fl (p, s) describing the coupled system of LAMS and SAVs is written as a product: q4(P, s) = IC/n(s;P)cpN(P).

(1)

The function cpN(p)describes 1 modes of LAMS in the state N = (N, , Nz. . . . , N,). p = (PI, /x2,. . . >pi) is a set of LAM coordinates. The function qn( s; p) depends on the dynamical SAV coordinates s = ( sI+i, s~+~,. . . , SJ)with the LAM coordinates p treated as parameters. The quantum numbers n = ( nl+, , n1+2,. . . , nJ refer to SAV modes. Both the 1c/,and the (ONfunctions can be determined as follows. First, the total vibrational Hamiltonian must be derived, H(P, s) = H,(s: P) + &(P; s) + T,,(P, s),

(2)

ff,(s; P) = TAs; P) + V,,(P, s),

(3)

H,(P; s) = T,(P; s) + I”(P)

(4)

where

and 2T, = P:G,,P,,

2T, = P,Tc,,P,,

2Tps = P;G,,P,

+ P:G,,P,.

(5)

Pi (where i = s or p) is the matrix of the momenta operators conjugate to the ith coordinates, i.e., Pi = -ifid/dqi. The vibrational tensors G,,, G,,, and G,, = G,, depend, generally, on all coordinates p and s. They will be considered in detail in Section IV. In the adiabatic method, the SAV wavefunction is derived from

[HAS; P) - &b)l~n(s;

P) = 0.

(6)

INTERNAL

LARGE

AMPLITUDE

MOTION

HAMILTONIAN

139

Let us note that H,(s; p) depends on p as a parameter, because it does not include the differential operator P, . Therefore the eigenvalue E,(p) is a function of p and can be interpreted as a contribution to the LAM potential function generated by the SAVs. Now, the LAM wavefunction can be found by searching for a minimum of the functional EN = (+;;‘(p, s) IH(p, s) I+i(p, s))~,~, provided that &,(s;p)is normalized, i.e., (tin(s; ~)lJ/,,(s; P)), = 1. Th e minimum condition leads to the Schrodinger-like equation for (Pi, [T,(P) + V’(P) + E”(P) -

-6vIevb) = 0.

(7)

where 2-,(P) = (IC/n(s; P)I T,(P, s) + T,,(P, s)l#&,

s)>s.

(8)

The operator ‘tT, can be written explicitly (by taking into account the definitions of T, and T,,) as T,(P) = $p;fG;,,p,

+ A,,(P) + &(P)

+ {P:&(P))

3

(9)

where

(10)

G,, = (1C/nIGppb+h>s and the adiabatic corrections are given by

(11) (12) and B,, = -$(WW,s}IrC/n)s.

(13)

The braces { - - - > in Eqs. (9)-( 13) indicate that the Pi act only within them. The corrections A,, and A,, are, usually, very small because $,, vary slowly with p> so P,& is small. BP, depends on the form of G,,. When G,, is constant with respect to s, B,, is zero. Thus, in the adiabatic approximation, the main part of the crosskinetic energy term T,, does not contribute to the total vibrational energy. This is a source of a significant error, if the tensor G,, is not small. Hence, in order to assure a reasonable accuracy of the adiabatic method it is necessary to find coordinates p and s such that G,, is negligibly small. The solution to this problem will be presented in Section IV. III. DETERMINATION OF THE LAM EFFECTIVE POTENTIAL BASED ON THE ADIABATIC METHOD

FUNCTION.

In various approaches to the LAM-SAV problem (3, 19-30) the same common basic model is used. The following simplifying assumptions define this model: 1. Only one LAM is considered. 2. The tensor G depends only on the LAM coordinate p = pI . 3. G,,(p) is zero.

140

PYKA,

FOLTYNOWICZ,

AND

MAKAREWICZ

4. G,,(p) has a diagonal form since SAVs are described by normal 5. The potential function VJp, s) is a quadratic function of s: s vo,(P, s) = C i=2

The effects neglected as perturbations.

1

+ i b(p>sf

ii(P

in this model such as anharmonicity

1.

The eigenvalue

+ ki(p)s!]

(14)

of the SAVs can be treated

Let us apply the adiabatic theory to this model. The energy E,(p) mode energies ei( p) determined from ( $[-h2Gii(p)d2/dsj!

coordinates.

+ li(p)Si - ei(p)}p,(si;

is a sum of one-

p) = 0.

(15)

ei( p) is equal to ei(p) = hwi(p)(% + 1/2) - wi(p),

(16)

where We

=

[GdP)k(P)l”2

cp; is the wavefunction from the equilibrium. A,,(P)

wi(P) = C(p)/2k(p).

and

(17)

of the harmonic oscillator displaced by A;(p) = --li(p)/k;(p) The adiabatic correction A,, calculated for our model is = C ]a;(p)Ml

+ hwi(~)Q~)(ni

+

l/2)1,

(18)

where ai = G,,[ hdln(wi/Gii)/dp]*/

16,

di = G,,( dAi/dp)2/2Gi;

(19)

and M, = n; + nj + 1. Now, the effective potential

for LAM can be written

?nn(P) = vo(P) - w(P)

+ C Ihiji(P)Cni

in the form + 1/2) + MiG(P)l,

(20)

where G(P) = w(P)[~

+

dib)l

(21)

and W(P) = C%(P).

(22)

Since vn(p) depends on a set of SAV quantum numbers n, its shape is different in various excited states of the SAV modes. We avail ourselves of this fact to propose an efficient method of determination of the potential parameters, from the experimental data. As an example let us consider a case when some of the SAV modes, say two modes s2 and s3, can be excited and transitions between LAM energy levels can be observed in these SAV states. Then, the LAM energy levels provide information about the

INTERNAL

LARGE

AMPLITUDE

MOTION

HAMILTONIAN

141

adiabatic LAM potentials v”‘,,,,(p) (the remaining nj = 0) in each excited SAV state, for example,

P,o = lJ” - W + )i(3& + &)/2

+ (3a2 + a3) + v:

PO, = I” - W + h(ij, + 3&)/2

+ (a2 + 3~) + l’:,

(23)

where v: = 2 (h&;/2 + a;) r>3

is a contribution from all remaining SAVs ( i > 3 ) . The one-dimensional vnZn,( p ) can be determined by fitting the series of the transition energies corresponding to various LAM levels in the states (n21/3). Thus, solving Eqs. (23) with. known vn,,,,( p), one obtains PO(P) = 2&o(P) - (l?o(P)

+ Qo,(P))/2.

h;,(P)

+

2azhJ)

=

Cob)

-

Gob)~

h&3(~)

+

2a3b)

=

do,

-

%dP)~

(24)

where VO(P) = VO(P) - W(P) + v:(P). As we can see, it is impossible to determine the pure LAM potential V’(p) due to the contributions of V:(p) and W(p) to v’(p). From the determined functions ;j,( p) + 2ai( p), the original functions ki( p) characterizing the molecular force field can be derived. The proposed fitting procedure is very simple and much more efficient than the traditional method based on solving the problem of several coupled vibrations and fitting the multidimensional potential function. The advantages of our method can be summarized in three points: (i) Only one-dimensional effective potential functions pn(p) are fitted, so onedimensional Schrodinger equations must be solved. (ii) The functions p,J p) are defined with a smaller number of parameters than the multidimensional potential V’(p) + V(p, s), so the fitting procedure is easier to perform. (iii) In each one-dimensional fitting only a subset of the experimental transition energies corresponding to a given excited SAV state is used. The above procedure allows for a partitioning of the whole fitting procedure into subprocedures which need much less computer time than the traditional approach. Let us illustrate how our method works. For testing purposes we have chosen the disilacyclobutane (DSCB) molecule, for which Killough et al. ( 14) determined a threedimensional potential function describing the ring-puckering LAM and two SAVs, viz. ring deformation and rocking of the hydrogen bonds. Their model Hamiltonian

142

PYKA,

FOLTYNOWICZ,

AND

MAKAREWICZ

fulfills all five conditions defining the basic model, formulated above. The ring-puckering potential function V” has been assumed as P(p)

= c2p2 + c4p4

(25)

and Vps(P,

s2, s3)

= (b2

+

C12P2)d

+ (b3

+

C13P2)&

(26)

= 13(p) = 0 and ki(p) = 2(bj + clip2). In order to determine the effective potential functions foe, Fro, and for, the experimental transition frequencies of DSCB taken from Table VI of Ref. (14) have been used. However, these data allow one to obtain pn with accuracy to an arbitrary constant co(n), i.e.. ?n = Vt;” + co(n). IJ’~ functions given in the form It means that in this model lo

c-F(P)

= c2(n)p2

+

c4(n)p4

have been obtained for the SAV states n = (n2, n3) = (0, 0), (1, 0), and (0, 1). The standard deviations of these fits are almost the same as those of Ref. ( Z4), so they are not reported here. The constants co(n) can be obtained from the condition co(n) = woo(n) + Eo(O) - Eo(n),

(27)

where woo(n) is the frequency of the transition from the state (N, n) = (0, 0) to (0, n), and EN(n) is the Nth LAM energy level calculated for the potential V:(p). The condition (27) is explained in Fig. 1, where the case of n = (0, 0) and IZ= ( 1, 0) is illustrated. Knowing wm( 1, 0) = 379.9 and wm(O, 1) = 437.7 (Ref. (14)), and calculating &(O, 0) = -46.29, &( 1, 0) = -63.82, and &(O, 1) = -39.24 (all quantities are given in cm-‘), we obtain the constants co(n) defining the potentials vn presented in Table I. Now, according to Eq. (24). it is possible to determine P’(p) (in the model of Ref. (14) only two SAVs are taken into account, SOV:(p) = 0) and Gi( p) + 2ai(p). We have found that the adiabatic corrections a2 and a3 are so small that they can be

V 00

FIG. 1. This figure explains how to establish the value of the co constant (see Eq. (27 )) with help of the effective potential functions & and PI;,, and the experimental transition frequency ww ( 1, 0)

INTERNAL

LARGE AMPLITUDE

MOTION HAMILTONIAN

143

TABLE I The Effective Potentials V&J) = co + c2p2+ c4p4(ckAre in cm-’ mdek) and hw, (in cm-‘) Determined for the DSCB Molecule for the Model of Ref. (14)

n = (r12,n3) CB

i

vo

=2

10+x

(0.0)

0.0

-8317

197.5

(1.0)

397.4

-9144

192.5

(0.1)

430.7

-78 16

194.6

0.0

-6 155

201.5

c4

hiA2

397.4

-626.6

-5.059

hZ3

430.7

501.1

-2.767

omitted without any decrease in accuracy. The functions &i(p) and v’(p), given as quadratic polynomials, are presented in Table I. In order to compare the results of Ref. (14) with ours, let us note that in the model ofRef. (14) l,(p) = 0, SO Ai = 0 and Gi = w, (see Eqs. (19) and (21)). Thus, we can obtain directly ki(p) as ( Lji)2/ Gil. The functions kl( p) and kj( p) are presented in Fig. 2. It is important to note that in our method, the forms of k,(p) are not assumed but come as the result of calculations. Conversely, in the traditional method of Ref. ( 14) the forms of ki( p) had been postulated and they were taken as quadratic functions of p. As follows from Fig. 2 our ki(p) cannot be well approximated by parabolas in a wide range of p. However, in the range -0.2 < p < 0.2 A, excluding the maximum of k,( p) they can be fitted by parabolas: ki(p) = 2( 6, + ?lip2) which yields 6, = 22.16 x 103, c”,2= -10.67 X 104, 6, = 6.836 X 103, and El3 = 1.304 X 104. We can see that these values are very close to the values bz = 22.3 X 103, cl2 = -10.96 X 104, b3 = 6.823 X 103, and cl3 = 1.359 X lo4 obtained in Ref. ( 14) (for the units, see Table V in Ref. (14)).

FIG. 2. Plots of the k,(p) and k3( p) functions defining the model of Ref. ( 14) obtained with help of the method presented in this paper (continuous lines) and in Ref. ( 14) (dotted lines).

144

PYKA,

FOLTYNOWICZ,

AND

MAKAREWICZ

Also, our potential I”(p) is very similar to that of Ref. (14) (in this case W(p) = 0 because 1i = 0). The barrier height calculated for our V’(p) is 82.5 cm-’ and the potential minima occur at p = 20.1422 A, whereas the analogous values of Ref. ( 14) are 80.8 cm-’ and &O.1424 A, respectively. In summary, the three-dimensional potential defined by P”(p) and ki( p), obtained by us, is practically identical to the potential of Ref. ( I#), if the same model of coupled LAM and SAVs is assumed. This fact proves the validity of our approach. Now, let us focus our attention on the LAM-SAV model used above. The assumed potential V( p, s2, s3) is incomplete because low-order terms such as Z13ps3(the symmetry types of p and s3 are the same) and l12p2s2 (~2 is totally symmetric, see Ref. (14)) were excluded. These terms hindered the solution of the eigenvalue problem when the variational Ritz method was used. In our approach these terms do not cause any difficulties. However, calculations performed with all these terms included in the potential lead to the conclusion that, in general, the inverse spectral problem does not have a unique solution. To explain this closer, let us consider two different potentials: V( p, s) = V’(p) + VPs(p, s) with VP, defined by Eq. ( 14), and V’(p, s) = V”(p) + V&(p, s) with vbs defined by Eq. ( 14) in which 1i(p) are replaced by I:(p) and ki( p) are unchanged. Let us choose V”(p) in such a way that VfO(P)

-

W’(P)

=

VO(P)

-

W(P).

(28)

Then, if the small adiabatic corrections aj and dj can be neglected, the effective LAM potentials v,,(p) and f”,(p) will be identical. This means that the total potentials V( p, s) and V’(p, s) will give the same eigenvalues. Hence, it is impossible to determine separately Zi(p) and V’(p) in the course of the fitting procedure, if the Z;(p) are not too large (so the bi are small). Thus, the inverse spectral problem for the LAM-SAV model has no unique solution. The conclusion presented above was drawn on the basis of the adiabatic approach. However, it is general and is not affected by the calculation method. In order to prove this fact, let us consider the simplified two-dimensional models including only the coordinates p and s3. For the first model we take V’(p) given in Eq. (25) and 13(p) = 0, and for the second model V”(p) = V’(p) + I:(p)/2k3(p) and &(p) = Ii3p. For both models k,(p) = 2( b3 + c13p2). The transition energies calculated variationally for both models are presented in Table II. As we can see, the transitions hardly change with the parameter 1i3. Even for the value of II3 approaching about 10% of b3, these transitions change only by 0.1 cm-’ or less. Thus, if the experimental transition frequencies have an accuracy of 0.1 cm-‘, as in the case of DSCB, the models with different II3 (but not too large) are indistinguishable. So, the inverse spectral problem for the LAM-SAV model does not have a unique solution also when the variational method of calculations is used. Another important question concerning the LAM-SAV models is the calculation of the tensor G. In the above models the cross elements G,,, according to Ref. (14), have been omitted. As we already stated in Section II, the elements G,, should be negligible to assure a good accuracy of the adiabatic method. However, for the threedimensional model of DSCB considered here, the cross term G1.3g 0.033 ( 17) is not small (the diagonal elements: G ,,, g 0.008 and G3,3 r 0.403). The influence of G1,3

INTERNAL

LARGE AMPLITUDE

MOTION HAMILTONIAN

145

TABLE II The LAM Transition Energies G, (in cm -’ ) Calculated for the Two-Dimensional LAM-SAV Model ( n3 = 0) with Various Values of the Parameter I,, (in cm-‘/A rad) N

113 = 0

350

100

0

4.35

4.34

4.30

1

52.74

52.77

52.85

2

35.28

35.24

35.15

3

51.23

51.22

51.20

4

56.73

58.72

56.68

5

62.62

62.61

62.58

Note. The remaining parameters are taken from Table V of Ref. ( 14)

on the transition energies can be easily estimated by applying the Van Vleck theory which gives the modified LAM element G-i,, g G,,, - GT.JG3.3. G,,, differs from G,J by 30%, so G1,3 considerably perturbs the calculated transition energies and should not be omitted. Moreover, all elements Gij for the three-dimensional LAM-SAV model in Ref. ( 14) have been calculated, in fact, for one-dimensional models (compare Table V in Ref. (14) with Table XII in Ref. (17)). Such an approach seems to be inconsistent. So, the derivation of a proper tensor G should be analyzed carefully. IV. CURVILINEAR

COORDINATES

AND THEIR ORTHOGONALIZATION

In order to solve the pure vibrational problem for a polyatomic molecule, proper internal coordinates q = (p, s) must be defined and the kinetic energy T must be expressed in terms of generalized momenta p conjugate to q. The general method of derivation of T proposed by Eckart (44) is very convenient in numerical realization and it has been recently developed by Laane and co-workers ( 16. 17). A general algorithm of this method is the following: (i) Define explicitly the Cartesian atomic vectors rk (whose coordinates are defined with respect to some molecular axis system) as functions of q. (ii) For a given q compute rk( q) and all derivatives drk( q)/aqi. (iii) Using rk and drk/dqi compute the inertial tensor I, the Coriolis interaction matrix X, and the pure vibrational matrix Y, which define the kinetic matrix K:

K=

I [ XT

x Y

1

so that

2T = (oT, qT)K

w 0q

.

(29)

For the definitions of the matrices I, X, and Y, see Refs. (16, 17, 44). (iv) Invert K to obtain the matrix G which allows the expression of T in terms of generalized momenta: 2T = (JT, pT)G

J

0 P

*

146

PYKA.

FOLTYNOWICZ,

AND

MAKAREWICZ

The w and J are the angular velocity vector and angular momentum vector, respectively. In order to derive the quantum operator ‘i’ corresponding to T, J and p should be replaced by the corresponding operators and the pseudopotential U(q) should be added to V(q) (45). Steps (ii)must be repeated for each value of q. In this way G(q) can be computed on a discrete set of points and then can be approximated by appropriate analytical functions, e.g., polynomials and periodic functions. We employed the above algorithm to numerically compute G(q) for four-membered ring molecules like DSCB. The functions rk( q) fully determined by taking into account all internal coordinates are presented in the Appendix. In the previous works not all coordinates were considered; for example, all vibrations breaking the C, symmetry of a molecule ( 16) or all nonsymmetric deformations of the ring (46) were neglected. Moreover, to calculate the effective masses for LAMS, only the LAM coordinates p were treated as dynamical and all SAVs were frozen out (s = 0). In such an approach, instead of the whole matrix K, only its submatrix

X,h 0) y&J, 0) 1

(31)

was formed and then inverted to obtain GR(p, 0) = [KR(p, O)]-‘. GR defines the rigid model of LAMS. However, another model of LAMS, the so-called static model, with frozen SAVs is possible (4 7). It is defined simply by Gs = G (p, 0). In this model the whole matrix K should be inverted and in the matrix G obtained the SAV coordinates should be fixed. These models are, in general, inequivalent. The rigid model depends on the dimension (i.e., on the number of LAM coordinates included in the model), whereas the static model does not. Let us compare the calculated matrices GR and GS for the DSCB molecule. In Fig. 3 the elements GF, and G?,, , referring to the ring-puckering coordinate p = ql, are

si=

0

0.15

1

O.Ul -1.0

G 1.1 -0.5

0.0 C

0.5

1.0

trod1

FIG. 3. Elements GfI, and Gs,I, referring to the RP coordinate p = q, plotted as functions of p. Numbers assignedtocurvesdefinetherigidmodels: 1,2. 3,and5 mean I-d(l).2-d(l,9). 3-d(l.9. 13),and5-d( I. 9, 13. 25, 29). respectively (the abbreviation “n-d( 1. 2, , WI)” means n-dimensional model including the (~1, pi. . , pm)coordinates). The symbol S refers to the static model.

INTERNAL

LARGE AMPLITUDE

MOTION HAMILTONIAN

147

plotted as functions of p. We can observe very significant differences between the rigid models of various dimensions. In all these models only the coordinates of symmetry B,, are included. In models where coordinates q, of different symmetry (in the D2h group) are also included, Gt, does not change considerably. It can be easily understood, because all nondiagonal elements GE are quite small (for q = 0 they are exactly zero). In the five-dimensional R model including all B,, coordinates, Gr, is very close to G?,i for small p. However, the difference between G & and G?,, is essential for large p. The explanation of this fact is simple. For p =# 0 the symmetry of the molecule becomes Cz,, for which four additional vibrations of the same type as p( A, ) occur. They couple with p and influence the form of GR. It is worth studying the dependence of GR and GS on isotopic substitutions. The results collected in Table III show that this dependence is different for the static and rigid models. For example, Gs,, is independent of the masses of atoms attached to the ring, whereas G& depends on them. For the first three R models GF, calculated at q = 0 is the same (for q # 0 there are small differences) because the coordinate p( B, U) does not couple with A, coordinates. Let us note that the matrix GS exhibits local properties. Its elements which refer to the ring coordinates qi-& are insensitive to isotopic substitutions of the atoms attached to the ring. Such locality is broken in the rigid models. Now, it is clear that the rigid models give different effective mass for a given LAM pi depending on whether the remaining LAM coordinates are taken into account or not. Such a dependence results from the fact that in the kinetic matrix K and, as a consequence, in the matrix G there are elements which couple a given coordinate p, with the remaining ones. Let us verify whether the elements Gri (i # 1 ), which couple the ring-puckering coordinate with the others, are small. At pl = 0 and s = 0 all elements Grj forj # 9, 13, 25, 29 vanish because the coordinates qj have symmetry types different than y(p)

TABLE III The Values G,(Z)/ G,,( N) Calculated at q = 0 for the normal (N ) and Isotopic ( I ) Modifications of DSCB Specified in the First Column The model

Static

Rigid l-d(l)

Z-d(l,Z)

0.1566 0.9710 0.9936

0.7566 0.9710 0.9938

3-d(1.2,3) 3-d(1.2.9) 30-d(l-30)

Gl l(I)/G1,l(N) DSCB-ki* DSCFJ-Cl3 DSCB-Sizs

0.7566 0.9710 0.9936

0.6747 0.9666 0.9903

1.0000 0.9461 0.9696

G2 2(I)/G2,2(N) DSCB-HP

0.9362

0.9344

0.9362

1.0000

DSCB-Cl3

0.9761

0.9332

0.9761

0.9226

0.9765

1.0000

0.9765

1.0000

DSCB-Si2'

1

Nofe. The abbreviations are the same as those in the legend to Fig. 3. DSCB-H 2 means that in the DSCB molecule all nuclei H ’ are replaced by Hz.

148

PYKA,

FOLTYNOWICZ,

AND

MAKAREWICZ

= B,, (see Table VIII). The remaining four elements for the coordinates of Bi, symmetry are illustrated in Fig. 4. As we can see, they are not small. For s Z 0 all elements Gij can be nonzero; however, as we checked, even for si = 0.3 they are very small. In order to eliminate the cross elements Gps, the coordinates q = (p, s) can be transformed to q’ = (p’, s’) so that G’,I,~ r 0. Since the tensor G can be written as (G);, = Cm;’ k

(32)

(aq;/&)(%Jark),

it transforms, under the substitution qi + q:, to G’ = AGAT,

(33)

Aik = dq:/dqk.

(34)

where

Now, let us write the coordinate transformation

in the following form: (35a)

P’ = P + C B&)x k St

=

s.

(35b)

Then, calculating A from Eq. (34) and using Eq. (33) we obtain GbP = aG,,aT

+ BGFsaT + aG,,BT + BG,BT

(36a)

G’, = G,

(36b)

GbS = aG,, + BG,,

(36~)

where a=l+c;

(c),j

=

c

(37)

(dBik/aPj)Sk*

i:l-_rI:--_i -1.0

-0.5

FIG. 4. Elements G,j (j refers to coordinates group) plotted as functions of p.

0.0

Q

belonging

[rod3

0.5

I.0

to the B,. representation

of the L&, symmetry

INTERNAL

LARGE

AMPLITUDE

MOTION

HAMILTONIAN

149

From the condition Gbs = 0 we can determine B: B = -aG,,G;‘.

(38)

However, this matrix depends on s (because a is a function of s) , so it does not have an assumed form B(p) . To avoid this difficulty, we can determine B for s = 0. Then a(p, 0) = 1 and B(P) = -G,,(P,

O)G,‘(p,

0).

(39)

Gbs calculated with such a matrix B(p) is zero for arbitrary p but only for s = 0. The nonzero contribution to Gbs is proportional to s and can be neglected in zeroth-order approximation. For s = 0 the matrix CL, simplifies to G;,(P’, 0) = Gpp(~, 0) - G&J, O)G,‘(p,

O)G%,

0).

(40)

Using the relation GK = 1, it can be proved that G;,(P’, 0) = ]KR(p,

W-‘.

(41)

This relation implies that G z, ( p’, 0 ) has the same functional form as G R( p , 0 ) defining the rigid model. However, the submatrix G:,(p’, 0) has exactly the same form as G,,( p, 0) (see Eq. (36b)) which is obtained by inverting the whole matrix K. The transformation B(p) can be written in terms of the matrix D = [ G,( p , 0)]-’ , where G, is the pure vibrational submatrix of G. Using the relation G,D = 1, we can prove that B(P) = I&,%.

(42)

For a one-dimensional LAM model p represents only one LAM coordinate and then D,-b = 1/D,,. In this case we obtain the transformation derived by Quade (30). He obtained his R transformation for a very special molecular axis system defined by nonholonomic conditions which make all Coriolis terms vanish. In a later work of Guan and Quade (48) the additional transformation ( T transformation) of molecular axes has been applied to eliminate the Coriolis terms originating from the SAVs. By performing both R and T transformations they obtained a block diagonal matrix K (in the zeroth-order approximation) which after inverting gives a G matrix with Gzs = 0. Our derivation requires no preliminary conditions to be put on the molecular axis system and is valid for an arbitrary number of LAM coordinates. Also, if computations are performed numerically, the form of B given in Eq. (39) is more convenient than that of Eq. (42)) because only two matrix inversions are needed to obtain G from K and G;’ from G,. According to Eq. (42), three inversions are necessary because the additional inversion must be performed to obtain D,-b from D,,. Naturally, for one LAM coordinate both forms of B are numerically almost equivalent. The SAV coordinates s’ can also be orthogonalized to eliminate the cross elements in the G, submatrix. The transformation PC= PI

(43a)

s# = LT(p’)s’

(43b)

150

PYKA,

FOLTYNOWICZ,

AND

MAKAREWICZ

leads to a new submatrix G&, which for s” = 0 takes the form Gzs = LTGksL.

(44)

L can be determined numerically as a matrix diagonalizing G&. Similarly, the matrix diagonalizing simultaneously G:, and the quadratic (with respect to s’) part of the potential P’bscan be found if I’bs is known. In order to orthogonalize the LAM and SAV internal coordinates, a computer program which operates in two stages has been written. In the first stage the matrix G(p, 0) is computed at several points pi (only one LAM coordinate is considered), according to the algorithm described above. Then, G,( p, 0) is inverted and the matrix B(p) is formed; see Eq. (39). Afterward, the elements B,~(P) are interpolated by polynomials of p. In the second stage of computations, the matrix G’( q’) is computed for a set of new coordinates q’ = (p’, 0); see Eqs. (35a) and (35b). Since a direct analytic dependence of the vectors rk on q’ is unknown, we must use a subroutine defining the vectors rk as functions of the original coordinates q. Thus, the q have to be computed for each value of q’ on the input. To find q corresponding to q’ the nonlinear equation (see Eq. (35a)) P = P’ -

c Blk(P)Sk k

is solved by using a simple iteration procedure which is rapidly convergent. In this way, all vectors rk and their derivatives with respect to the q: can be computed to form the matrix K’( q’). As a result, G’( q’) is obtained by inverting K’( 4’). The correctness of the above numerical procedure has been tested by computing G:@(p) from the rigid model. The elements GFP(p) and Gb,(p’) calculated for the same points p = p’ differ by less than 1OP4%,which proves the validity of our numerical method. The program computing G’ can be slightly modified to also orthogonalize the SAV coordinates. Only a standard subroutine diagonalizing G& must be added. We also propose another, simpler method allowing for the orthogonalization of all internal coordinates, by performing only one transformation: q’ = Liq.

(45)

Lo is chosen in such a way that the transformed matrix GL = L;fG(qO)LO is diagonal at the point qb = L;fqO. Also in the vicinity of this point, the nondiagonal elements of G L are small, if the matrix G(q) varies slowly with q. The transformation LOTgives local orthogonal coordinates q’ in the sense that GL is exactly diagonal at qb and almost diagonal in a quite large domain around qb. We tested and found that the Gtj calculated for the DSCB molecule are very small in a wide range of p and even of s also. In Fig. 5, some diagonal elements of G” and GL calculated for the DSCB molecule are compared. We see that they are almost equal in a large range of p. Thus, the simpler local transformation Lo of all coordinates can be applied instead of the two transformations B(p) and L(p). Naturally, this is not true when the matrix G(q) strongly depends on q .

INTERNAL

LARGE AMPLITUDE

MOTION HAMILTONIAN

0.0

Cj’ [rod1

0.5

151

1.0

FIG. 5. Plots of some diagonal elements G” (continuous lines) and G’- (dashed lines) calculated for the DSCB molecule.

V. RESULTS FOR DISILACYCLOBUTANE

Here, two LAM-SAV models of the DSCB molecule are analyzed. In the threedimensional model, the LAM ring-puckering (p) and SAV ring deformation ( s2) and in-phase rocking ( sg) are taken into account. In the four-dimensional model the SAV out-of-phase rocking (sia) is also considered. However, in order to derive a proper tensor G” expressed in orthogonal coordinates q”, all internal coordinates have to be included in the kinetic energy operator. The elements G$ are approximated by polynomials in the LAM coordinate, Gyi(p”) = i

gi( k)p””

(46)

k=O

whose coefficients a(k) are presented in Table IV. Using G;,, (p”), the effective potentials v,,( p”) are determined for the ground and some excited states of the SAVs observed in the infrared spectra ( 14). The results of one-dimensional fits are collected in Table V and the vn(p”) obtained are given in Table VI.

TABLE IV Diagonal Elements G;,(p”) (in amu-’ rad-‘) Calculated for the Orthogonal Coordinates q” as Functions of the Transformed Ring-Puckering Coordinate p”; see Eq. (46) Gii

g(0)

S(2)

E(4)

0.00232

a(6)

G1.l

0.07022

-0.01603

G2.2

0.21307

0.00065

0.00060

G9,9

0.54335

0.01751 -9.01114

0.00266

G10 10

0.52779

0.01466 -0.00196

0.00016

0.00030 -0.00064

152

PYKA, FOLTYNOWICZ,

AND MAKAREWICZ

TABLE V The Ring-Puckering Transition Energies GN (in cm-‘) of the DSCB Molecule in Various States n = ( n2, ng, n,,,) of the SAVs

state N

(1.0.0)

(0.0,0) exp.

exp. -cal 0.0

exp.

exp. 4.1

0.0

53.2

0.0

56.5

24.9

0.7

34.6

0.2

34.4

1.6 2.2

2

31.6

3

49.5

4

54.5

-0.8

5

61.0

-0

6

66. 0

-0.4

64.1

0. 4

67. 0

7

70.6

-0.2

10.2

2.0

8

74.7

9

78.5

-0.8

exp.-Cal.

0.4

56.0

0.0

exp.

0.6

3.1

1

4

exp. -cal.

64.8

0

-0

(0,B.l)

(0.1.0)

exp. -cal.

1.9

0 9

50.5

-0.4

53.0

55.5

-0.7

55.0

62. 0

-0.2

61. 1

-1.0

0. 0

68. 0

-1

71 5

0. 1

70 7

-0.7

0. 0

75.4

0.1

75.5

0.2

79.5

0.7

78.8

48.6 51.5 4

57.4

-2.4 0.3 -1.4

-1.0

0

0.3 -0.2

10

82.0

0.4

82.5

0.2

11

65.3

0.6

85.5

0. 1

12

86.2

0.6

13

91.2

0.9

14

93.9

1.1

, Note. The experimental values are taken from Ref. (14).

Now, the potential function p’(p”) and the force constants k;(p”) can be determined by applying the adiabatic theory presented in Section III. In the case of the threedimensional model, the barrier height of @ii is 82.9 cm-’ (compare with 80.8 cm-’ for Laane’s model of Ref. (14)) and the potential minima are at p” = 0.4205 rad which corresponds to 0.14 16 A for Laane’s coordinate (compare with 0.1424 A for Laane’s model). The analogous potential p’pv of the four-dimensional model has a barrier height 84.2 cm-’ at 0.4229 rad (0.1423 A), respectively. The potentials @ii and pf: differ by hw,o/2, so they have slightly different barriers. The polynomials approximating v” and who; (i = 2,9, 10) for the four-dimensional model are presented in Table VI. The force constants k,(p”) are plotted as functions of p” in Fig. 6. The obtained potential I/(p”, s’;, s $, s’;~) can be expressed in terms of the original coordinates p and Si, as V = V’(p)

+ i 7 [fT(p)si + ky(p)s? I

+ C kgsis,] + higher order terms. J

New kp can be expressed by parameters of V(p”, s”) and coefficients of the transformation q --* q”. The diagonal force constant kq is plotted as a function of p in Fig. 6 also (the remaining k? differ slightly from ki and are not presented in the figure). This constant is compared with kz calculated with Gz,* determined from a one-di-

INTERNAL

LARGE

AMPLITUDE

MOTION

153

HAMILTONIAN

TABLE VI The Effective Potentials v,Jp”) = c0 + c2p”* + c4pw4 (ckAre Given in cm-’ and hw, (in cm-‘) Determined for DSCB ” = ol*.“s,nla)

i

n

ca

c2

C4

(0.0.0)

0.0

-959.2

2627

(1.0.0)

397.5

-1057.5

2601

(0,l.O)

430.7

-903.7

(0.0.1)

436.4

Go “O2 Gg Y0

0. 0

2601

-952.3

2669

-941.3

2631

397.5

-98.33

-22.26

430.7

55.46

-26.50

436.4

Note. P” and hwk are determined

6.942

for a four-dimensional

rad-k)

41.41

model.

mensional rigid model, as in Ref. ( 14). We can see that k2 differs considerably from kq, which means that rigid models for the SAVs are inadequate. VI. SUMMARY

In this work a fully numerical approach to the internal LAM-SAV Hamiltonian of polyatomic molecules is developed. It is applied to four-membered ring molecules in which all internal motions are described by curvilinear coordinates. The kinetic energy tensor G(q) is constructed numerically. The rigid models of the ring-puckering motion are considered and the tensor G R( p ) is calculated for them. It is shown that GR( p) depends on a number of the internal

I.?

C [rod1

FIG. 6. Plots of the k*(p”), kg(p”), k,,,(p”) force constants (continuous lines) calculated for the DSCB molecule with the help of the four-dimensional model. The dotted line represents function k,(p) calculated with the use of Gs2 obtained from the one-dimensional rigid model. For graphical presentation purposes, k2 ( p ) is multiplied by 0.1. The dashed line is a plot of the k;( p ) force constant function (see the text )

1.54

PYKA,

FOLTYNOWICZ,

AND

MAKAREWICZ

motions which are taken into account in the model. This results from the fact that LAM and SAV coordinates are, in general, nonorthogonal. The transformation of LAM and SAV coordinates to new orthogonal coordinates q’ is derived (see Eqs. (35) and (39)). For these coordinates the nondiagonal block G,,(p’, 0) = 0; i.e., all LAM-SAV zeroth-order cross terms in the kinetic energy vanish. Also, all zeroth-order cross terms between SAV coordinates can be removed. An alternative transformation which gives local orthogonal LAM and SAV coordinates for which the tensor G is diagonal at a given point is considered. It is proved that this transformation can also be practically useful. To separate SAVs and LAMS, the adiabatic theory is employed and the effective LAM Hamiltonian is obtained. It has a simple form if the SAVs are described by harmonic oscillators. Their effective masses and force constants can depend on LAM coordinates. The corrections to the LAM Hamiltonian due to anharmonicity can be easily calculated analytically using the perturbation theory. In Table VII, our approach is compared with the HBJ and Quade theories (only the pure internal motions are considered ) . A simple and efficient method of derivation of the SAV-LAM potential function from the experimental transition energies is proposed. It is shown, however, that the inverse spectral problem has no unique solution. The three- and four-dimensional potential functions of DSCB molecule are derived and compared with the potential function of Ref. ( 14). It is shown that one-dimensional rigid models used to derive the effective masses for SAVs give, in general, incorrect force constants. APPENDIX

So far, a model allowing for the treatment of an arbitrary number of LAMS in ring molecules has not been investigated. Here, we will consider all internal motions of the DSCB molecule for which a three-dimensional model has been recently considered (14). It is a typical four-membered ring molecule and other molecules of this kind can be considered as modified DSCB with some atoms replaced by other atoms. The numbering of the hydrogen (H) and ring atoms of DSCB is shown in Fig. 7. TABLE VII Comparison

of the LAM-SAV

Theories

(3), Quade (JO), and This Work

Quade

HBJ 1.

of HBJ

This

work

SAV coordinates

2. SAV-LAM

Rectilinear

Curvilinear

Curvilinear

Sayvetz

R transformation

B transformation

Van Vleck

Adiabatic

ortho-

gonalization

condition

3. SAV-LAM separation 4. SAY

Van

Vleck

theory

LAM

dependent

One

LAM

theory

theory

basis

functions 5. Number

of LAMS

Nofe. HO, harmonic

oscillator.

HO

LAM

independent

One

LAM

HO

LAM dependent Several

LAMS

HO

INTERNAL

LARGE AMPLITUDE

MOTION HAMILTONIAN 9

5 1

1 6

3 L

12

v

155

7

2

10

0

FIG. 7. Numbering of the atoms of the DSCB molecule. The Si atoms have the numbers 1 and 3. and C atoms have the numbers 2 and 4. H atoms have the remaining numbers,

The internal coordinates describing the ring modes, which can be treated as LAMS, are illustrated in Fig. 8. They are chosen as follows: R,are the lengths of the vectors Ri, R is the length of the ring diagonal, and p is the angle between two planes shown in the figure. This angle describes the ring-puckering motion and it is assumed to be positive (negative) if Y, is positive (negative). This coordinate is different from the traditionally used distance between the two ring diagonals. However, we prefer p since the atomic coordinates are simple functions of p. To describe internal motions of H atoms, we use the global system of axes (GSA) and local system of axes (LSA) connected with the ring atoms, similarly as in ( 16. 46). The GSA and one of the LSA are illustrated in Fig. 8. The GSA is defined as follows: (i) the Z axis lies in the plane spanned by the second, third, and fourth ring atoms and is directed to the third atom.

FIG. 8. The GSA (X, Y, Z) and LSA (xr, yr, z,. ) placed at atom I. The polar-spherical coordinates (d, 19, cp) of the fifth H atom are shown in the upper frame. The instantaneous configuration of the four-atomic ring is described by the internal coordinates p, R; and R. y, and yz are auxiliary quantities.

PYKA,

156

FOLTYNOWICZ,

AND

MAKAREWICZ

(ii) the X axis passes through the second and fourth atoms and is directed to the fourth atom; (iii) the Y axis is perpendicular to the X and Z axes and is chosen in such a way that the (X, Y, Z) GAS is right-handed. Each LSA is based on three ring atoms. Let us consider the LSA placed on the first ring atom. The xL and zL axes lie in the plane formed by two vectors: RI and R4. The zL axis bisects the angle fi between the vectors RI and -R4, so it is parallel to the unit vector e3 =

(W&

The yL axis is parallel

- RI/&)/~

cos(P/2),

IR;l.

to the unit vector e2 = R4 X R,fR4R,sin

and so is perpendicular vector

Ri s

where

p,

to the (RI, R4) plane. The xL axis is determined

by the unit

el =e2Xe3. Now, the position of each H atom can be determined by its local vector rk (in Fig. 8 the r: is shown) with LSA coordinates (x,, yi, zi). However, they are inconvenient as internal coordinates and instead of them we use the polar spherical coordinates (d,, tY;, pi) which are related to (xi, yi, Zi) by the following equations: xi = disin y,

=

TACOS

tpi

disk dish (pi

zi = dices 29i. The position of the ith H atom, measured from the origin of the GSA, is rc = r” + RF, where RjL is the vector of the jth ring atom which is the center of an LSA and r” = el,xi + eyyi + e3jzi. Let the GSA coordinates of the basis vectors ekj be (XE, YE, Z;)j, then the GSA coordinates (X7, Y?, Zp) of rf; can be expressed as

In order to express (X0, Y?, Zp) in terms of internal coordinates, we have to use the above equations and write the vectors ek, and RF as functions of the ring internal coordinates. The GSA coordinates of the ring atoms as functions of Ri, R, and p are given by XF = D, - D2,

Yf = D4sin p

Xi

Z: = -D4cos

= -D2,

Xk = R - D2,

Z;=D3

p

INTERNAL

LARGE AMPLITUDE

MOTION HAMILTONIAN

157

TABLE VIII Symmetrized Internal Coordinates Describing Vibrational Modes of the DSCB Molecule Symmetry %?I

Vibrational mode

=2" Ring vibrations P AR

92

AR1

93

q

AR2

+

94

RF' Rdef

Rmod

95

AR3

96 AR4 SiH Wi) 97

A85

scis "es

.

96 99 410

Ae9 *B10

rock

CH Wi) 911

A@7

scis he6

*

912

A811

913

A812

914

rock SIH @vi) @5

915

% @9

.b q16 917

twist

wag

@lo 916 CH (api) &7 @8

919 920

twist

921 922

wag

4 @ll %2

SiH (di> % Ad6

923 *

Ad9

924

str

925

Ad10 926 CH (di) M7 M6 %I Ml2

927 .b

'28

str

929 q30

Note. Abbreviations used: RP, ring puckering; Rdef, ring deformation: Rmod, ring modes; scis, scissoring; rock, rocking; etc. ip (op), in (out of) phase; s, symmetric: a, antisymmetric; and * means the transformation A.

158

PYKA, FOLTYNOWICZ, AND MAKAREWICZ

and the remaining coordinates are zero. The D; are defined by D, = R,cos

Dz

yI,

=

D4 = Risin yI,

03 = Rzsin y2,

RZCOS72,

where yI = arc cos[(R2 + R: - Rz)/2R,R

1

y2 = arc cos[(R’

I.

and + R: - R:)/2RzR

Now, the GSA components of the vectors R; can be deduced immediately, for example, R, = R: - Rk = (-D, , - D4sin p, D&OS p). These components are required in order to calculate the angles p,j between Ri and -Rj and to determine the vectors ekj. Therefore, the atomic GSA coordinates can be expressed by the internal variables qL = (p. R, {R;},

{di, oi, CP~}).

However, the origin of the GSA does not lie at the center of mass of the molecule, so we must shift all the atomic GSA coordinates by (XCM, YCM, ZCM) where, for example, XCM = 2 m,XySAf 2 m;;

XPSA = Xf- or XC .

The above introduced internal variables qL have a local character. They are inconvenient when the molecule has a high symmetry. However, they can be transformed to collective symmetrized coordinates q. The transformation qL,to q can be described by an orthogonal matrix A of the form

which transforms a chosen set of four local coordinates qL( 4) = to four symmetrized coordinates q( 4), for example,

[I [I 93 94 95 96

MI

=A

2;

;

ARi = Ri - (Ri),.

AR4

The inverse transformation is given by qL( 4) = Aq( 4)) because A-’ = A. The point groups D2,, and C,, have been previously applied (19) to classify the vibrational modes of DSCB. In the solid (gas) state, the DSCB molecules are planar (puckered), so the D2,, ( C2,) point group is correct. For our purposes both symmetry groups will be useful. The forms of the symmetrized coordinates of DSCB and their symmetry types with respect to the DZh and C2” groups are presented in Table VIII. All the introduced internal coordinates are curvilinear and they can adequately describe LAMS. RECEIVED: February

19, 1990

INTERNAL

LARGE

AMPLITUDE

MOTION

HAMILTONIAN

159

REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. Il. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48.

C. C. LIN AND J. D. SWALEN, Rev. Mod. Phys. 31,841-892 (1959). W. R. THORSON AND I. NAKAGAWA, J. Chem. Phys. 33,994-1004 ( 1960). J. T. HOUGEN, P. R. BUNKER, AND J. W. C. JOHNS, J. Mol. Specfrosc. 34, 136-172 ( 1970). P. R. BUNKER AND J. M. R. STONE, J. Mol. Spectrosc. 41, 310-332 ( 1972). P. R. BUNKER AND B. M. LANDSBERG,J. Mol. Spectrosc. 67,374-385 ( 1977). CH . V. S. R. RAO, J. Mol. Spectrosc. 89, 197-2 13 ( 198 1). D. PAPOUSEK,J. M. R. STONE, AND V. SPIRKO, J. Mos. Spectrosc. 48, 17-37 ( 1973). D. C. MOULE AND CH. V. S. R. RAO, J. Mol. Spectrosc. 45, 120-141 (1973). M. KREGLEWSKI, J. Mol. Struck 60, 105-108 (1980). P. JENSENAND P. R. BUNKER, J. Mol. Spectrosc. 94, 114-125 (1982). L. A. CARREIRA, R. C. LORD, AND T. B. MALLOY, JR., in “Topics in Current Chemistry” (M. J. Dewar et al., Eds.), Vol. 82, pp. l-95, Springer-Verlag, Heidelberg, 1979. M. A. HARTHC~CK AND J. LAANE, J. Phys. Chem. 89,4231-4240 ( 1985). L. E. BAUMAN, J. M. COOKE, J. R. VILLARREAL,P. M. KILLOUGH, AND J. LAANE, J. Phys. Chem. 86, 2000-2006 (1982). P. M. KILLOUGH, R. M. IRWIN, AND J. LAANE, J. Chem. Phys. 76,3890-3898 ( 1982). P. M. KILLOUGH AND J. LAANE, J. Chem. Phys. 80,5475-5480 ( 1984). J. LAANE, M. A. HARTHCOCK, P. M. KILLOUGH, L. E. BAUMAN, AND J. M. COOKE, J. Mol. Spectrosc. 91,286-299 (1982). M. A. HARTHC~CK AND J. LAANE, J. Mol. Spectrosc. 91, 300-324 ( 1982). R. M. IRWIN AND J. LAANE, J. Phys. Chem. 82,2845-2850 ( 1978). A. R. HOY AND P. R. BUNKER, J. Mol. Spectrosc. 52,439-456 ( 1974); 74, l-8 (1979). P. JENSEN,P. R. BUNKER, AND A. R. HOY, J. Chem. Phys. 77,5370-5372 ( 1982). P. JENSENAND P. R. BUNKER, J. Mol. Spectrosc. 99, 348-356 ( 1983). P. R. BUNKER, Annu. Rev. Phys. Chem. 34,59-75 ( 1983). P. JENSENAND P. R. BUNKER, J. Mol. Spectrosc. 118, 18-39 (1986). R. BEARDSWORTH, P. R. BUNKER, P. JENSEN, AND W. P. KRAEMER, J. Mol. Spectrosc. 118, 40-49 (1986). P. JENSEN, Comput. Phys. Rep. 1, l-56 (1983). V. SPIRKO, J. M. R. STONE, AND PAPOUSEK,J. Mol. Spectrosc. 60, 159-178 ( 1976). V. SPIRKO, Mol. Phys. 38, 1761-1766 (1979); 43, 723-735 (1981). V. SPIRKO, J. Mol. Spectrosc. 101, 30-47 ( 1983). A. SAYVETZ, J. Chem. Phys. 7,383-389 ( 1939). C. R. QUADE, J. Chem. Phys. 65, 700-705 (1976); 73,2107-2114 (1980). D. H. CRESS AND R. C. QUADE, J. Chem. Phys. 67,5695-5701 ( 1977). B. I. STEPANOV, Nature (London) 157,808-814 (1946). A. WITKOWSKI, J. Chem. Phys. 47,3645-3648 ( 1967). J. BRICKMAN, Int. J. Quantum Chem. 7,427-436 ( 1973). P. RUSSEGERAND J. BRICKMANN, J. Chem. Phys. 62, 1086-1093 ( 1975). N. D. SOKOLOV AND V. A. SAVE’EV, Chem. Phys. 22,383-395 ( 1977). J. MAKAREWICZ AND A. WIERZBICKI, Chem. Phys. Lett. 108, 155-160 ( 1984). B. R. JOHNSON, R. T. SKODJE, AND W. P. REINHARDT, Chem. Phys. Lett. 112,396-402 ( 1984). H. R. DUBAL AND M. QUACK, J. Chem. Phys. 81,3779-3791 ( 1984). J. MAKAREWICZ, J. Mol. Spectrosc. 105, l-7 ( 1985). J. MAKAREWICZ, A. WIERZBICKI, AND J. KOPUT, Chem. Phys. 97,31 l-319 ( 1985). B. R. JOHNSONAND W. P. REINHARDT, J. Chem. Phys. 85,4538-4556 ( 1986). T. CARRINGTON, JR., J. Chem. Phys. 86,2207-2223 ( 1987). C. F. ECKART, Phys. Rev. 47, 552-568 ( 1935). R. MEYER AND H. H. GONTHARD, J. Chem. Phys. 49, 15 lo- 1520 ( 1968 1. V. SZALAY, J. Mol. Spectrosc. 110, 172-173 (1985). J. MAKAREWICZ, J. Mol. Spectrosc. 130,316-336 (1988). Y. GUAN AND C. R. QUADE, J. Chem. Phys. 84,5624-5638 ( 1986).