8922
J . Phys. Chem. 1993,97, 8922-8928
An Efficient Method for Optimal Modes Analysis of Vibration: Application to HCCD Karen E. Aiani and John S. Hutchinson’ Department of Chemistry and Rice Quantum Institute, Rice University, Houston, Texas 77251- 1892 Received: February 10, 1993; In Final Form: May 24, I993
The quantitative analysis of vibrational spectra of polyatomic molecules in the high-energy regime requires a determination of the modes which optimally describe the vibrational motions of the nuclei at these energies. Observed vibrational spectra in small polyatomics indicate substantial regularity in the vibrational motion, implying that such a set of “optimal modes” should exist, but the experiments do not provide a direct means of characterizing these modes. We present here a computationally efficient theoretical method for optimal modes analysis of exact multimode vibrational eigenstates. This algorithm consists of direct numerical integration of selected projection coefficients which reveals the extent of zeroth-order character of these eigenstates. Application of this method is presented for the analysis of high-energy vibrations in monodeuterioacetylene.
I. Introduction
have been presented for performing such an analysis, principally including adiabatic basis ~ a l c u l a t i o n s , coordinate-dependent 2~~~~ Vibrations of highly excited polyatomic molecules play sigSCF calculations,2~*9and natural orbital calculations.30Jl nificant roles during chemical reactions, but a quantitative Our own approach, as presented in previous paper~’~-3~ and description of the relationship of detailed vibrational dynamics reviewed in section 11, is to gauge the quality of a particular to reactivity remains underdeveloped. This circumstance should coordinate description by calculating the zeroth-order character be expected to change rapidly, however, due to the development of exact multidimensional eigenstates. Briefly, a set of zerothof increasingly sophisticated experimental methods of observing order states is uniquely defined for each coordinatesystem chosen. the photophysical and photochemical properties of vibrationally In our analysis, the coordinates which produce the “best” basis excited molecules. These new methods include stimulated are defined to be the optimal modes. More specifically, in the emission pumping,14 vibrationally mediated photodissociation,5 optimal coordinates which describe a given exact eigenstate, one vibrational predissociation ~pectroscopy,6~~ double-resonance member of that set of zeroth-order states has the largest possible s p e c t r o ~ c o p y ,and ~ ~ ~ultrafast laser photochemistry.14.15 In projection onto the eigenstate. In our earlier work, on cycombination with increasingly well-resolved single-photon abanoacetylene22hydrogen cyanide,33and acetylene,” we calculated sorption methods,’”’* these experiments provide significantly these projections directly and exactly by diagonalization of the more data about the nature of vibrationally excited states than have been traditionally available. Hamiltonian matrix in a basis of zeroth-order states. One of the most remarkable aspects of the new data is that the This optimal mode analysis is similar in spirit to determining observed spectra reveal significant assignable structures, which the optimal modes based on a natural orbital expansion, and a is evidence of approximate regularity in the high-energy vibrations discussion of the relationship between these approaches is of small molecules. Thus, the data can be interpreted via a set presented in section V. We note here, however, that the optimal of zeroth-order states in terms of which the high-energy spectrum modes we found in our study of HCN are the same modes can be assigned. As a consequence, it is possible to infer from determined by using the natural expansion method. In fact, we the spectra the extent of mixing between zeroth-order state~~~320 have found that for HCN the natural orbitals are insignificantly and, from this, the strength of the energetic coupling between different from the one-dimensional eigenfunctions.33 these states. In polyatomic molecules with significant densities Analysisof the optimal modes of vibration by any of the methods of states, the intramolecular relaxation lifetimes of zeroth-order listed above requires a knowledge of the essentially exact states can be deduced;21*22 in vibrational predissociation spectra, vibrational eigenstates,preferablyvia diagonalizationor recursion one can infer zeroth-order state-specific unimolecular reaction of a large Hamiltonian, for each choice of coordinates; furtherrates.61 more, optimization of the coordinates mayrequirea large number On the other hand, what cannot be inferred directly from the of choices. Thus, although previouswork has illustrated theutility spectral data is the nature of the modes which defines these zerothof an optimal mode analysis via an expansion technique, these order states. Specifically, the quality of the zeroth-order methods are clearly impractical for routine analysis of molecules assignmentsin the high-energyregime reveals that an approximate with more than three vibrational degrees of freedom. separation of vibrational coordinates must exist. The coordinates In this paper, we present an extension of our previous approach, in which this separation is most accurate form the set of optimal incorporating a new computational algorithm35 for performing modes of vibration. The experimental spectra, however, do not the optimal mode analysis. The new algorithm, which is based reveal these coordinates. It is one thing to accurately assign a on a highly accurate approximation derived from the HEG vibrational band as, for example, 3vl + 2v2; it is quite another method,3638 requires only a single diagonalization of the to determine the motions of the atoms which define v 1 and v2. In vibrational Hamiltonian and is thus sufficiently fast as to permit the sense that the more recent double-resonance studies of highroutine analysis for molecules with many degrees of freedom. In energy vibrational states contain information about the extent of mixing betweenvibrationalmodes?+Itit is all the more significant section 11,we present a developmentof our new method of optimal that one determine the definition of these modes. mode analysis. The method is rigorously tested for accuracy by application to hydrogen cyanide in section 111, for which an exact Accurate quantitative interpretations of these excited vibrational motions require determination of the molecule’s natural, optimal mode analysis has already been performed. In section or “optimal”, modes. The need for such an optimal mode analysis IV we present an extension of this technique and apply it to a has been recognized in the literature, and a variety of approaches three degree of freedom molecule, monodeuterioacetylene. 0022-3654/93/2097-8922%04.o0/0
0 1993 American Chemical Society
The Journal of Physical Chemistry, Vol. 97, No. 35, 1993 8923
Efficient Method for Optimal Modes Analysis 11. Development of the Method
required integral explicitly by projecting a single basis state onto
A. Criterion for Optimal Modes. Consider for illustration a two-mode coupled oscillator system. (The results are straightforwardly extended to more than two modes.) In terms of some set of coordinates S1 and S 2 , we can write the Hamiltonian as
*/:
H(S1,SJ = Hi(S1)
+ Ht(S2) + H,2(S,,S2)
(1)
UlUZ
Of course, by orthonormality of the zeroth-order states, (w~l(wllul)lu2)= bwlulbwz, eq 5 becomes
(~2I(WllJI/) = ~(W,,W,),/ where HP is simply that part of the Hamiltonian which depends only on mode j . The definition and choice of S1 and S 2 are arbitrary; we have assumed nothing at this point about the proper zeroth-order nature of the Hamiltonian. The zeroth-order basis defined by this choice of coordinates is Iul)lu2) where
This basis can be used in a linear variational calculation of the two-mode eigenstates, $1,
H(S,,S,) */(SlS2) = W / ( S l S 2 )
(3)
by expanding the multimode eigenstates over the zeroth-order states
(4) UI
02
If the choice of coordinates S1,S2actually separated the twomode Hamiltonian, then we would conclude that these are the optimal coordinates for the molecularvibrations. Quantitatively, this is indicated when, for each choice of I , there is a single coefficient, C(u,srz),/, which is equal to one. More realistically, the optimal coordinates will not exactly separate the Hamiltonian, but might nearly do so such that a single coefficient will be close to one for each choice of 1. We can thus quantitatively measure the “quality” of a coordinate description by calculating the coefficients, C(ul,u2),/, and looking to see how near to one the dominant coefficient is for each choice of 1. We can also examine these coefficients as we vary the choice of coordinates to find the set which best describes the Hamiltonian. We have previously employed this analysisto determine the optimal vibrational modes of several small polyatomic molecules.32-34 It is relevant to note here that the zeroth-order expansion described above is not defined in same manner as the natural orbital expansion method, although there are apparent similarities. This distinction was alluded to in the introduction and will be expanded upon further in the discussion section. B. Direct Integration of C(a,ri)~. The most straightforward means to determine the required coefficient C(ul,u2),~ is via linear variational theory. Specifically, for each choice of the set of coordinateswe can find the zeroth-order vibrational Hamiltonian. For each of the single degrees of freedom in that Hamiltonian, we can find the zeroth-order eigenfunctions by solving the onedimensional Schrbdinger equation. We can then calculate the multimode Hamiltonian matrix in a basis of zeroth-order product functions, (u1)Ju2). Diagonalization of this matrix then provides directly all Ncoefficients for all N eigenvectors. This approach, used in our previous studies, is clearly accurate, but it is highly inefficient and therefore excessively time-consuming. Most importantly, for a single eigenstate, we are interested in only one or a few coefficients, C(u,,uz),/, for a few zeroth-order states. By contrast, the diagonalization of H produces NZ coefficients,of which all but a few are disregarded entirely in the analysis. Thus, it is much more efficient, although less straightforward, to determine only those coefficients required in the analysis. This can be done by direct integration. We can write the
(6)
A single coefficient, C(wl,wd,/,is thus a single multidimensional integral. Our method for performing this integration is based on our extensionn of the transformation method, introduced by Harris, Engerholm, and Gwinn3”3* (the HEG method). The HEG method has been shown to be mathematically equivalent to Gaussian integration,37 although our extension of it yields an expression for solving problems of a different form. For purposes of the optimal mode analysis, it is instructive to write explicitly the desired integral as
c(wI,wz),/ = (W2(~2)l(Wl(~l)lrl/(~1,~2))
(7)
Note that we have written the two-dimensional wave functions explicitly in terms of two sets of coordinate pairs, (Rl,R2) and (SI,&). (R1,Rz) represents a possible set of optimal modes for the two-dimensional vibrations. The goal of the analysis is to choose (R1,R2)so that C(wl,w,),/ is as close to 1 as possible; hence, they are chosen to separate the Hamiltonian as much as possible. By contrast, the set (SI,&) can be chosen strictly for computational convenience, rather than physical significance. Any choice of coordinates can be made here. We prefer to use internal modes (bond-angle coordinates) for (S132) since the vibrational Hamiltonian is most simply written in these coordinates. Note also that the full two-dimensional eigenfunctions,$1, are only needed in the coordinates (Sl,S2), and hence the multidimensional problem needs only to be solved once in whatever coordinates are most convenient. is thus determined as in eq 4. Only the zeroth-order eigenfunctions are required in the coordinates (R1,R2). 1 . Calculation of Iul(Sl))Ju2(S2))and J.1. To determine the zeroth-order eigenfunctions, Iu,(Sj)) as in eq 2 (j= 1,2), we use the HEG method. First, we define a primitive basis of harmonic oscillators {&(SI))in the coordinate S,,and we use these to determine an auxiliary basis {xk(S,))which later facilitates the numerical integration. The x k basis is defined by calculating the matrixof theoperator Sjin the harmonicoscillator representation and then diagonalizing it. Specifically, we calculate
Xnm = ( $nISj#m ) and diagonalize the resultant matrix
(8)
A = T‘XT to define the new basis functions X k
(9)
We now adopt the set of functions { x k j for the linear variational basis for determination of Iu,(S,)). The HEG approximation for the potential energy matrix elements gives The kinetic energy matrix elements are calculated by transformation from the harmonic oscillator representation:
After combining the potential and kinetic energies to construct the one-dimensional Hamiltonian matrix and, subsequently, diagonalizing it, we can very efficiently determine the zeroth-
8924 The Journal of Physical Chemistry, Vol. 97, No. 35, 1993
order eigenfunctions as superpositions of the functions {xk)
Aiani and Hutchinson can rewrite the target integral
C(wl,wl),l in
eq 7 as
(w2(R2)1(w1 ( ~ l ) l $ / ( S l 3 2) )=
Having found the zeroth-order eigenfunctions, we now use them as a basis for determining the multi-mode eigenfunctions, +/(Sl,S2).Two items are of note here. First, calculation of the Hamiltonian matrix elements in the Iu1)lu~)basis is greatly facilitated by having each zeroth-order eigenfunction expressed as a superposition of the functions { X k ] , since we can implement the HEG approximation in two dimensions. Second, for the same reason, use of the ( X k ] basis simplifies and expedites calculation of the target integral, C(u1,9),1. This will be demonstrated below. By combining eqs 4 and 13, we thus have the double sum
where
The integral inside the summation = ~ ~ ~ , ( R ~ ) ~ ~ ~ ~ l ( R ~ ) ~ X k l ~ S(23) ~))~x~l(~
can be viewed as an integral over either (SI&)or (Rl,Rz).The two sets of coordinates are, of course, functions of one another. We take the approach that the integral is over (SI&),and R1 and R2 are evaluated as functions of SIand Sz. The twodimensional integral here is thus not separable into two onedimensional integrals. However, evaluation of this integral is greatly facilitated, as mentioned above, by the presence of the functions X k in the integrand. We have recently shown3s that a one-dimensional overlap integral containing X k may be very accurately approximated by m
The ak, ,,values are determined when finding the one-dimensional zeroth-order functions as in eq 13. The C~uI,ul),~ values are determined when diagonalizing the full two-mode Hamiltonian in the zeroth-order basis defined by the (S&) coordinates, as in eq 4. This completes the right half of the integrand in eq 7. 2. Calculation of (w2(Rz)l(wl(R1)J. The left half of the integrand in eq 7 consists of one element of the zeroth-order basis in the "new" coordinates (Rl,R2). This can be determined in a manner exactly analogous to that described above for the (SI,&) coordinates. The set of functions thus determined is different, of course, since the zeroth-order partitioning of the two-mode Hamiltonian is different in the new coordinates. For purposes of clarity of definition of required quantities, we repeat the development of the previous section, this time, for the (R1,Rz) coordinates. Just as before, we begin by defining a set of primitive harmonic oscillators yn(R,) for each coordinate R j These are then used to determine an auxiliary basis ?)k(Rj),which, as we mentioned earlier, make solution of the target integral much simpler. As before, we first calculate
>
Ynm ( YnlRjlYm and diagonalize the resultant matrix
B = U'YU
(16)
(17)
This defines the new basis functions 7lk
In this equation, 8 ( x ) is any arbitrary function satisfying the same boundary conditions as the primitive harmonic oscillator basis &o(x). The factor TOk is one of the coefficientsdefining the function X k in eq 10, and the value Akk is the eigenvalue of the matrix X in eq 9 corresponding to the eigenstate X k . Hence, an overlap integral involving X k requires only two functional evaluations. This is easily extended to two (or more) degrees of freedom by performing the integral over S I first and then over Sz, The result is that the required integral M in eq 23 can be approximated very efficiently and accurately by
where we have assigned S1 = Aklkl and S2 = AkDkl. Combining this result with eq 22 completes the desired integration of the optimal mode criterion, C(wl,wl),l. We are left only to show that, in actual applications, this approach is both sufficientlyfast and sufficiently accurate for routine optimal mode analysis of vibration. This is demonstrated in the next section for our test molecule, hydrogen cyanide, for which we have previously determined the optimal mode coordinates. Section IV then presents an extension of the method to HCCD, which has three degrees of vibration.
n
Constructing the one-dimensional Hamiltonian as before and diagonalizing, we find that
fori = 1,2. Substituting in the definition of 9 from eq 18 yields an expression for the zeroth-order functions in (R1,Rz)in terms of the harmonic oscillator primitive basis:
where
This completes the left half of the integrand in eq 7. 3. Calculation O ~ C ( ~ ~By , ~combining ~ ) , J . eqs 20 and 14, we
IXI. Demonstration of the Method: Application to HCN As our first application of the direct integration method presented above, we consider the stretching vibrationsof hydrogen cyanide. The HCN molecule has often served as one of the standard prototypes for spectroscopicobservationand theoretical calculation of high-energy vibrational spectra and eigenstates. Smith et 01.39 observed the high overtone spectrum with very good resolution. The stimulated emission pumping spectrum of vibrationally excited HCN was observed by Yang et aL49@ Accurate quantum calculations of HCN vibrational eigenstates have been presented by BaEiC et al.41.42and by Brunet et ~ 1 . ~ Our own workon HCN, as previouslyreported>3demonstrated the necessity of interpreting the overtone spectrum in terms of an energy-dependent set of optimal modes. As such, this is an appropriate system for testing the method just described. As in our earlier study, we make use of a very accurate potential surface constructed by Botschwina" which predicts the high-energy band
3
The Journal of Physical Chemistry, Vol. 97, No. 35, 1993 8925
Efficient Method for Optimal Modes Analysis
TABLE I: Comparison of HCN Eigenstate Projection Coefficients Calculated by Diagonalization and by Integration state
matrix diagonalization bond modes optimal
direct integration optimal
(091 ) (092) (093) (0~4) (03)
0.964 0.933 0.904 0.878 0.854
0.996 0.993 0.989 0.984 0.977
0.993 0.989 0.984 0.979 0.972
(190) (290) (3,O)
0.964 0.922 0.870
0.993 0.982 0.963
0.990 0.979 0.960
(192) (291) (193) (292) (194)
0.762 0.729 0.679 0.568 0.606
0.962 0.944 0.943 0.905 0.922
0.960 0.942 0.942 0.905 0.922
positionsvery precisely. His surface describes only the stretching vibrations and is based on ab initio electronic structure calculations in thecoupled electron-pair approximation (CEPA-1). Recently, Carter et al.45 and Gazdy and Bowman46 have used variational calculations tooptimize the potentialof HCN tocreateanaccurate three-mode global potential for the molecule. However, for our model incorporating only the two stretching modes, we have used Botschwina's surface without modification, as it reproducs the observed single photon absorption spectrum sufficiently well. Our earlier study determined that coordinate rotation of the bond mode axes produced a set of optimal coordinates in which the Hamiltonian was nearly separable. This was achieved by taking linear combinationsof bond mode coordinatesand applying the criterion described in section 11. For each different choice of coordinates, the full multidimensionalHamiltonian expression was diagonalized to find the complete set of eigenstate projection coefficients. As stated in a prior section, rediagonalizing the problem for each new set of coordinates is very accurate but is also very time-consuming and inefficient. In this case, we have done so only to verify the accuracy of our multidimensional integration method. Toapply thenew integrationmethod within themodelof HCN, we first solve the full multi-dimensionalproblem by diagonalizing the system as expressed in some "simple" set of coordinates, (271,s~). We chose this to be an "unrotateds (bond mode) coordinate set. Then, using the rotations determined from our previous study, we solve only for the one-dimensional zerothorder eigenfunctions in coordinates (&,&). Finally, we utilize the integration scheme outlined above to find the corresponding zeroth-order projection coefficientsfor the particular eigenstates of interest. Table I compares the values of the projection coefficients calculated by these two methods: complete matrix diagonalization and direct numerical integration. Both techniques use a primitive basis of 25 harmonic oscillators for each mode and a multidimensional product function basis with 10 functions from each mode. The coefficients calculated by our integration approach are in excellent agreement with those found by the traditional linear variational methods. The numerically integrated values are never in error by more than 1% (assuming the diagonalized values to be exact), and the errors areusually considerablysmaller than that. Table I also demonstrates the utility of the optimal mode description of the HCN vibrations, but we refer the reader to our previous paper on the subject for a full discussion of these results.33 Even more importantly, the time required to perform the integration by the functional evaluation method is trivial, particularly when compared to the computationally intensive matrix diagonalization technique. For this two-mode system, a single projection integral requires only about 5 s on a SPARC-
TABLE Ik Comparison of HCCD Eigenstate Pro'ection Coefficients Calculated by Diagonalization and by fntg(ntioa matrix diagonalization optimal
bond modes
state
direct integration optimal
0.872 0.144 0.617
0.980 0.952 0.913
0.979 0.944 0.907
0.969 0.934 0.887 0.819 0.728 0.668
0.999 0.994 0.977 0.947 0.916 0.903
0.999 0.993 0.975 0.941 0.914 0.901
0.856 0.738 0.637 0.552
0.990 0.982 0.97 1 0.959
0.995 0.984 0.972 0.957
0.257 0.202 0.653
0.913 0.870 0.948
0.9 10 0.863 0.950
+,
station 1 while a matrix diagonalization requires about 15 min on the same machine. Having verified that our technique is both sufficiently fast and accurate for calculation of projection coefficients, we now apply this method to the problem of a larger molecule for which no feasible scheme was previously available for analyzing its optimal modes.
IV. Extension of the Method A. Application to HCCD. As a more extensive test of our integration technique, we consider the stretching vibrations of monodeuterioacetylene. Due to its linear geometry and fouratom structure, HCCD is a logical successor to the hydrogen cyanide molecule. In addition, there have been a variety of spectroscopicand theoretical studies performed on acetylene and various isotopically substituted acetylenes. Halonen et uL4' have presented a very accurate model for acetylene which accurately reproduces the known overtone band origins. BerryI9 has presented an analysis of the intensities and band positions of the high-lying overtones of acetylene. More recently, Lehmann" presented a detailed study of the intensities and Smith et al.49 examined the overtones of a bendlstretch system. Also, the stimulated emission pumping spectrum of vibrationally excited HCCH was observed by Sundberg et a1.2 We have chosen to use the monodeuterated molecule, HCCD, to avoid complications arising from the symmetry of HCCH. For this three-mode system, our target integral is
In a method analogous to that in two dimensions, this leads to a summation over the following integrals
= ( Yn3(R3)1(7b(R2)1 ( Y n , ( R 1 ) l X k , ( S 1 )
) k k , ( S 2 ) )Ixk3(s&
)
(27) Finally, we apply our accurate approximationfor overlapintegrals of this type to yield the following expression
The above direct integration technique is compared to the full diagonalizationmethod in Table 11. We findvery good agreement of the projection coefficientscalculated by the integration method with the exact (linear variational) values. Errors are all less than 1% in magnitude, with the average being 0.4%. What is striking
0926
Aiani and Hutchinson
The Journal of Physical Chemistry, Vol. 97, No. 35, 1993
TABLE III: Optimization of Selected States for HCCD ~
bond modes eigenstate (3,0,0) (0,2,1)
state 13)10)10) 12)11)10) 10)11)12) 10)13)10) 10)10)13) 10)12)11) 11)12)10)
~
~~~~
optimal modes
coefficient
state
coefficient
0.887 0.403
13)10)10) 12)11)10)
0.969 0.052
0.603 0.568 0.377 0.257 0.224
10)12))1) 10)13)11) l0)ll)ll)
0.916 0.276 0.230
in this molecule is the dramatic superiority of the rotated coordinate description over the bond mode coordinates. Clearly, the rotated coordinates do a much better job of separating the vibrational Hamiltonian, and thus of approximating the optimal coordinates, than do the bond modes. As with HCN, 25 harmonic oscillators were used in the primitive basis, and a multidimensional product function basis was constructed with 10 functions from each mode. With three modes, the reduction of the computational time when using the direct integration method for optimal mode analysis is even more dramatic than it was for the HCN system. Direct integration of a single coefficient only requires about 10 s (on a SPARCstation l+), whereas the full matrix diagonalization takes over 3 h on the same machine. This affirms our belief that the direct integration approach provides an effective and efficient means for determining the optimal modes of a molecule. We can expand on the results in Table I1 to demonstrate the utility of an optimal mode analysis of high-energy vibrational states. Table I11 shows the results of an optimal mode analysis of two selected eigenstates. For the eigenstate (3,0,0), the zerothorder state 13)10)10)in the bond modes has a fairly large projection coefficient of 0.887. However, the contribution from zero-order state 12)11)10), which has a coefficient of 0.403, is still fairly substantial. Using our new approach to calculate projection coefficients, we found a set of optimal coordinates for which the state 13)10)10) has coefficient 0.969. An even more obvious indicator of the improvement, however, is the large decrease in the 12)11)10) coefficient to 0.052. Thus, in the optimal mode coordinates,the eigenstate (3,0,0) reveals itself to be only slightly mixed in character. Indeed, it appears that, in the optimal coordinates, the state 13)10)10) is very nearly an eigenstate of H, indicating that the Hamiltonian nearly separates in these coordinates for this energy distribution. The eigenstateproperly labeled as (0,2,1) appears greatly mixed in character when analyzed in terms of bond modes, as shown in Table 111. No single zeroth-order state clearly dominates, and the zeroth-order state with the correct nodal arrangement, )0)12)11),is well down the list of contributing states. Upon rotating to the optimal coordinates, however, we see that the zeroth-order state 10)12)11) describes the eigenstate very well. In the bond modes, there is no way to properly assign this eigenstate, given the small coefficient, 0.257, of zeroth-order state 10)12)11). Rotating to a more optimal coordinate set, though, allows for simple assignment of the state as (0,2,1), due to its dominant coefficient of 0.9 16. The reduction of contributions from other zeroth-order states also implies that we have a greatly improved description of this eigenstate, and, correspondingly, of the vibrational Hamiltonian. B. Description of Optimal Mode Character. We have shown in the above sections that we can indeed find a set of optimal modes for both HCN and HCCD which more nearly separate the Hamiltonian than the normal modes. However, if the determined optimal modes were only minimally different from the normal modes, there would be no real significance to the findings. In order to determine how the found optimal mode coordinates vary with increasing energy, we have plotted essentially “fourdimensional”contour plots for each of the eigenstates of the H-C overtones.
Figure 1. “Contour” plot of the wave function for an eigenstateof HCCD at 18 603 cm-l assigned as (6,0,0). Lobes composed of dark dots are of positive sign while light colored dots are negative.
Y
W Y Lo
u
~~
CH stretch
Figure 2. Contour plot of (6,0,0) eigenstate of HCCD presented through one face of the cube so that only the mixing between the modes of CH stretch and CC stretch is displayed.
These plots consist of a three-dimensional axis system (for the three stretching local modes) with the different values of the wave function represented by different colored dots. Dots of one color, then, make up one energy “contour” in three dimensions. Because of the three axes, the overall plot ends up looking like a cube, with different colored “blobs” defining the maxima of the wave function. Nodes of the wave function occur where regions of different colors are adjacent to one another. We can choose to look at the cube through just one of its faces and determine the extent of mixing, demonstrated by the amount the plotted wave function is rotated relative to the perpendicular bond axes. Figure 1 contains a plot of the unrotated wave function for an eigenstate of HCCD at 18 603 cm-l assigned as (6,0,0). From Table 11, we see that the 16)10)10) zeroth-order state has a coefficient of just 0.668 (44% character). By examination of the plotted wave function, however, it is obvious that the eigenstate is well-defined, but not in terms of the chosen axis system. In Figure 2 the plot is presented so that only one face, representing the C-H stretch and C-C stretch, is displayed. If there were no coupling between these modes, the plot should lie parallel to one or both of the axes. Upon viewing Figure 2, however, this is clearly not the case. The contours appear to be “pivoted” about the middle so that they now lie at an angle to the axes. By examining the wave function plotted in this manner, we can determine the extent to which these particular modes are coupled together. Examination of the cube from other faces allows us to measure the extent of mixing between other pairs of the modes.
Efficient Method for Optimal Modes Analysis
l5 10
14 .-.-. 5000
The Journal of Physical Chemistry, Vol. 97, No. 35, 1993 8927
I
-.-.-.
-*
10000
15000
20000
Energy(cm-1)
1
20
5000
10000
15000
20000
Energy(cm-1)
Figure 3. Angle of rotation away from bond mode axes for C-H overtone eigenstatesof HCCD for (a) C-H as mixed with CC stretch and (b) C-H as mixed with CD stretch. Dashed line in both figures corresponds to the associated normal-mode rotation angle. By using this information, we have been able to find a more optimal coordinate system for which the projection coefficient of this state then improves to 0.903 (82% character). The degree to which the axes are rotated is illustrated more completely in Figure 3. In this figure, the angles by which the individual axes are rotated (relative to the bond mode axes) are shown for the first seven excitations of the C-H stretch. Clearly, the axesvary in a nonsystematicmanner. Evenmoreimportantly, these optimal rotations are significantlydifferent from the normalmode rotations, even at low energies. From these observations, it is obvious that an optimal mode analysis is imperative for complete understanding of the mode character at each energy. V. Conclusions and Comments In this paper, we have presented what we believe is the first optimal mode algorithm which is sufficientlyaccurate and efficient for routine spectral analysis of molecules with severalvibrational degrees of freedom. Such an analysis of high-energy vibrational spectra in terms of an optimal set of vibrational modes has been pursued previously in the literature. We wish to comment here on the relationship of our approach to previous approaches, in terms of both the efficiency of the algorithms and their underlying formal principles. Our method for analysis of vibrational eigenstates, namely an examination of the zeroth-order expansion as a function of coordinates, is not the same as the methods presented elsewhere. However, it is most similar to the natural orbital (or natural expansion) meth~d.~O.~l The latter method is based on an analysis of the exact eigenfunction expanded over an optimally chosen set of product basis functions. There are two points to clarify in comparingour approach to the natural expansion approach. First, the expansions analyzed are not the same, since the zeroth-order basis is defined differently than the natural orbital basis. Second, the natural orbitals represent the optimal product basis only in a given set of coordinates, since they are not canonically invariant with the choice of coordinates.50.51 As such, in both approaches, we are examining a particular expansionof an exact eigenfunction, and that expansion depends on the choice of coordinates. The choice which optimizes either expansion is the set of optimal
modes. As mentioned in the introduction, we have shown in previous work33that the optimal coordinates and states calculated with both approaches are indistinguishable. We thus compare these approaches in terms of algorithmic efficiency. Prior to this work, the previous methods employed to optimize vibrational modes, including our own work and natural orbital calculations, necessitated rediagonalization of a multidimensional Hamiltonian for every choice of coordinate systems. Our new method requires only one vibrational Hamiltonian diagonalization,making the optimal mode search procedure much more feasible for larger molecules. We expect that the accuracy and efficiency of the new algorithm, combined with the previously demonstrated utility of optimal mode analysis, will make this technique very useful as a routine analytical tool for vibrational analysis. Furthermore, the method presented in section I1 is actually fairly simple to apply. It does assume that, in one coordinate system, it is possible to find essentially exact high-energy vibrational eigenstates. Although we have presented the eigenstate calculation in terms of zeroth-order states in the unrotated coordinates, this is not required. Any set of coordinates chosen for efficiency of calculation may be used. All that is required is an expansion of the multimode eigenstates over basis states which are eigenstates of the matrix of the position operator. This is the heart of the HEG method; it is furthermore the heart of the more recently described and popularly applied “discrete variable representation” (DVR).52 As such, the method of analysis presented aboveisreadilyinmrporatedintolargebasiscalculations using DVR, such as those presented by Wyatt et al.43 In fact, our method should prove valuable in the analysis of these highly accurate but quite complicated multidimensional eigenstates. As applied in this paper, our method of optimal mode analysis presumes that a set of optimal modes exists in the energy range near the eigenstate of interest. Of course, for chaotic vibrational motion at high energy, which is particularly likely in larger polyatomic molecules, a set of coordinates will not exist such that the vibrational Hamiltonian is even approximately separable. However, we expect that, even in this circumstance, an optimal mode analysiscould provevaluable. Recent work in our groups3~s4 and elsewhere5s*56indicatesthat, even when thevibrational motions of polyatomic molecules are chaotic, there can exist localized high-amplitude motions which persist for comparatively long periods of time. These motions generally encompass a limited set of the vibrational modes of the polyatomic, indicating an approximate separability which is useful on a time scale of many vibrational periods.5’ Furthermore, these localized motions can be expected to dominate the high-energy spectrum.5s As such, it is likely, in our view, that high-energy spectra, even in a quasicontinuum, can be fruitfully analyzed via optimal mode analysis. An additional possible application of this technique might be within the field of electronicspectroscopy. In analyzingelectronic spectra, it is often necessary to take into account the Duschinsky effect,Sg in which vibrational normal coordinates of different electronic states are displaced and/or rotated with respect to one another. This feature requires solution of an overlap integral of harmonic oscillators similar to that presented and solved in this paper. Kupka and Cribbm presented a very detailed means of solving these integrals by using lengthy generating functions. Alternatively, Chen and Pei61 applied an expansion method based on addition theorems of harmonic oscillator wave functions and associated Laguerre polynomials. We expect, however, that our direct integration method proposed aboveshould find use in solving problems of this nature, as it is a straightforward and efficient way to calculate this type of overlap integral. One additional benefit of the method proposed in this paper is that it permits calculations in nontrivial coordinate systems. The correct set of optimal modes describing a high-energy vibration may be (and perhaps should be expected to be) fairly
8928 The Journal of Physical Chemistry, Vol. 97, No. 35, 1993 complicated curvilinear coordinates. For example, Colbert and Sibert have optimized curvilinear stretching-bending coordinates in a model of C-H motion in methyl halides.62 Calculatingmatrix elements of the Hamiltonian expressed in such coordinates in a basis of zeroth-orderstates in these coordinatesproves to be both difficult and time-consuming. With our new method, exact eigenstates of the vibrational Hamiltonian need not be found in these more complicated coordinate sets: only the zeroth-order states are required, and these will be comparatively easy to calculate. Assuch, themethod presented here shouldbeof general utility for a wide variety of vibrational motions in polyatomic molecules.
Acknowledgment. This work was supported in part by grants from the Robert A. Welch Foundation of Houston, TX, and the National Science Foundation. Acknowledgment is made to the donors of the Petroleum Research Fund, administered by the American Chemical Society, for partial support of this research.
References and Notes (1) Abramson, E.; Field, R. W.; Imre, D.; Imnes, K. K.; Kinsey, J. L. J . Chem. Phys. 1985,83,453.
(2) Sundberg, R. L.; Abramson, E.; Kinsey, J. L.; Field, R. W. J. Chem. Phys. 1985,83,466. (3) Dai, H.-L., Ed.J . 0pt.Soc.Am.E 1990,7,1801containsa collection of recent results from SEP. (4) Yang, X.; Rogaski, C. A,; Wodtke, A. M. J . Chem. Phys. 1990,92, 2111. ( 5 ) Likar, M. D.; Baggott, J. E.; Sinha, A,; Ticich, T. M.; Wal, R. L. V.; Crim, F. F. J. Chem. Soc., Faraday Trans. 2 1988,84,1483. (6) Butler, L. J.; Ticich, T. M.; Likar, M.; Crim, F. F. J . Chem. Phys. 1986,85,2331. (7) McGinley, E. S.;Crim, F. F. J. Chem. Phys. 1986,85, 5741. (8) Luo, X.; Rizzo, T. R. J. Chem. Phys. 1990,93,8620. (9) Luo, X.; Fleming, P. R.; Seckel, T. A.; Rizzo, T. R. J. Chem. Phys. 1990,93,9194. (10) Luo, X.; Rizzo, T. R. J . Chem. Phys. 1991,94,889. (1 1) Fleming, P. R.; Li, M.; R i m , T. R. J. Chem. Phys. 1991,94,2425. (12) Fleming, P.R.; Li, M.; Riuo, T. R. J . Chem. Phys. 1991,95,865. (13) Fleming, P. R.; Rizzo, T. R. J. Chem. Phys. 1991.95, 1461. (14) Scherer, N.F.; Doany, F. E.; Zewail, A. H.; Perry, J. W. J . Chem. Phys. 1986.84, 1932. (15) Scherer, N. F.; &wail, A. H. J . Chem. Phys. 1987,87,97. (16) McIlroy, A,; Nesbitt, D. J. J . Chem. Phys. 1989,91,104. (17) McIlroy, A.; Nesbitt, D. J. J . Chem. Phys. 1990, 92,2229. (18) Kerstel, E. R. T.; Lehmann, K. K.; Mentel, T. F.; Pate, B. H.; Scoles, G. J. Phys. Chem. 1991, 95,8282. (19) Berry, M. J. Vibrational Photochemistry and Photophysics In Proceedings of the Robert A. Welch Foundation Conferences on Chemical Research; Welch Foundation: Nov 1984;Vol. XXVIII, p 133.
Aiani and Hutchinson (20) Dubal, H.-R.; Quack, M. J. Chem. Phys. 1984,81,3779. (21) Sage, M. L.; Jortner, J. Adu. Chem. Phys. 1981,47,293. (22) Heller, D. F.; Mukamel, S . J. Chem. Phys. 1979,70 (1). 463. (23) Johnson, B. R.; Reinhardt, W. P. J. Chem. Phys. 1986,85,4538. (24) BaZiC, Z. J. Chem. Phys. 1991,95,3456. (25) Gerber, R. B.; Ratner, M. A. Ado. Chem. Phys. 1988, 70,97. (26) Thompson, T. C.; Truhlar, D. G. J. Chem. Phys. 1982,77, 3031. (27) Lefebvre, R. In?. J. Quantum Chem. 1983,23, 543. (28) Moiseyev, N. Chem. Phys. Lett. 1983,98,233. (29) BaEiC, Z.;Gerber, R. B.; Ratner, M. A. J. Phys. Chem. 1986, 90, 3606. (30) Moiacyev, N.;Schatzverger, R.; Froelich, P.; Goscinski, 0.J. Chem. Phys. 1985,83,3924. (31) Holmer. B. K.; Certain. P. R. J. Phvs. Chem. 1985.. 89.. 4464. (32) Hutchinson, J:S. J . Chem. Phys. 1h5,82,22. (33) Fleming, P. R.; Hutchinson, J. S.J. Chem. Phys. 1989,90, 1735. (34) Fleming, P. R.;Hutchinson, J. S.Comput. Phys. Commun. 1988,51, 59. (35) Aiani, K.E.;Hutchinson, J. S.Comput. Phys. Commun. 1992,69, 46. Engerholm, G.G.;Gwinn, W. D.J. Chem.Phys. 1965, (36) Harris, D. 0.; 43, 1515. (37) Dickinson, A. S.;Certain, P. R. J. Chem. Phys. 1968,19,4209. 138) Shore. B. W. J. Chem. Phvs. 1973.59.6450. i39j Smith, A. M.; Jargensen, U. G.;Lehmann, K. K. J. Phys. Chem. 1987,87,5649. (40) Yang, X.; Wodtke, A. M. J. Chem. Phys. 1990,93,3723. (41) BaEiC, Z.;Light, J. C. J . Chem. Phys. 1987,86,3065. (42) MladenoviC. M.; BaEiC. Z. J. Chem. Phys. 1990.93,3039. (43) Brunet, J.-P.; Friesner, R. A.; Wyatt, R.E.; Leforestier, C. Chem. Phys. Lett. 1988, 153,425. (44)Botschwina, P. J. Chem. Soc., Faraday Trans. 2 1988.84, 1263. (45) Carter, S.;Handy, N. C.; Mills, I. M. Philos. Trans.R. Soc.London, Ser. A 1990,332,309. (46) Gazdy, B.; Bowman, J. M. J . Chem. Phys. 1991, 95,6309. (47) Halonen, L.; Child, M. S.; Carter, S . Mol. Phys. 1982,47, 1097. (48) Lehmann, K.K. J. Chem. Phys. 1989,91,2759. (49) Smith, B. C.; Winn, J. S . J . Chem. Phys. 1988,89,4638. (50) Moiseyev, N.; Wyatt, R. E. Chem. Phys. Lett. 1986,132,396. (51) Mokyev, N.; Friesner, R. A,; Wyatt, R. E. J . Chem. Phys. 1986, 85,331. (52) BaEiC, Z.;Light, J. C. Annu. Rev. Phys. Chem. 1989,40,469. (53) Hutchinson, J. S.;Marshall, K. T. J. Chem. Soc., Faraday Trans. 2 1988,84, 1535. (54) Marshall, K. T.;Hutchinson, J. S.J. Chem. Phys. 1991.95,3232. ( 5 5 ) Llorente, J. M. G.; Hahn, 0.; Taylor, H. S.1.Chem. Phys. 1990,92, 2762. (56) Holme, T. A.;Levine, R. D. J. Chem. Phys. 1988,89,3379. (57) Fumero, J.; Hutchinson, J. S.Unpublished research. (58) Levine, R. D.; Berry, R. S.J. Chem. Phys. 1989,90,2071. (59) Duschinsky, F. Acta Physicochim. URSS 1937,7,551. (60) Kupka, H.;Cribb, P. H. J. Chem. Phys. 1986,85,1303. (61) Chen, K.;Pei, C. Chem. Phys. Lett. 1990,165,523. (62) Colbert, D. T.; Sibert 111, E. L. J . Chem. Phys. 1989,91,350.
