4464
J . Phys. Chem. 1985, 89, 4464-4412
Acknowledgment. This work was supported in part by grants from the Scientific Affairs Division of NATO (RG. 138.81), the donors of the Petroleum Research Fund, administered by the American Chemical Society, and the National Science Foundation (CHE-8218216). We are grateful to Professor Michael Kasha
for bringing this problem to our attention, and to Dr. Aage E, Hansen for stimulating and helpful discussions. Registry No. I, 577-85-5; II(N), 13400-26-5; 111, 97877-63-9; 1V. 496-63-9; VI, 636-38-4.
Natural Orbital Analysis of Vibration-Rotation Wave Functions Bruce K. Holmert and Phillip R. Certain* Department of Chemistry, University of Wisconsin-Madison, Madison, Wisconsin 53706 (Received: April 8, 1985)
Natural orbital analysis has been applied to vibration-rotation wave functions of a model triatomic van der Waals system. Model potentials for both linear and T-shaped molecular geometries, along with various degrees of anisotropy, were considered. Calculations were not restricted to the ground state but included states with nonzero total angular momentum and excited bending vibrations. The vibration-rotation wave function was determined variationally by using a body-fixed Hamiltonian and then transformed and expressed in terms of natural orbitals. The natural orbital expansion of the wave function converged more quickly than the original expansion and thus facilitated the interpretation of the wave function. For the model potentials used, the natural orbital expansion consisted of one dominant term plus several correction terms. Through the use of the harmonic limit of the body-fixed wave equation, the nodal structure and spatial distributions of the correction terms were attributed to either vibration-rotation coupling or potential anharmonicity.
1. Introduction Advances in the spectroscopy of small and medium-sized have spurred the development of a variety of computational methods2v3which are used to calculate vibration-rotation spectra. One of the computational methods used for triatomic molecules employs the variational principle and assumes that the vibration-rotation wave function can be expressed as a sum of products of single variable vibrational functions and symmetric top rotational functions.* Alternatively, close-coupling methods have been applied to triatomic van der Waals molecules in which the set of close-coupled equations are solved by either explicit integration: perturbation t h e ~ r yor , ~basis set expansion.6 In all of these approaches, emphasis has been placed on the calculation of energies and other expectation values, with only minimal consideration of methods for analyzing the wave function. In this paper, the emphasis will be placed on the natural orbital analysis of variationally obtained wave functions in order to recover as simple a picture as possible of rotation-vibration states which are far from the harmonic oscillator-rigid rotor limit. An abundance of tools exist for the analysis of electronic wave functions: natural natural bond orbitals,I0 localized molecular orbitals,” etc. Similar tools for analyzing vibrationrotation wave functions have been introduced to build a greater understanding of molecular vibrations; however, application of the ideas of natural orbital analysis to nuclear motion has been limited. A model problem of two coupled harmonic oscillators was solved in closed form for its natural orbitals.12 Bishop and Cheung performed accurate nonadiabatic computations for H2+ and its isotopesi3 and then submitted the resulting wave functions, which involved both nuclear and electronic coordinates and consisted of over 300 terms, to natural orbital a n a l y ~ i s . ’ ~Stratt, Handy, and Miller used natural orbital analysis to examine the question of quantum chaos in a coupled oscillator system.15 Sage and Williams have used the analysis to study local modes in a coupled Morse oscillator system.16 None of these applications considered the rotational degrees of freedom. The purpose of the present contribution is to demonstrate the utility of natural orbital analysis when applied to vibration-rotation wave functions. Present address: Miller Fellow, Department of Chemistry, University of California-Berkeley, Berkeley, CA 94720.
0022-3654/85/2089-4464$0 1.50/0
To illustrate the application of natural orbital analysis to vibrational problems we studied a model triatomic van der Waals system and considered both linear and T-shaped molecular geometries with various degrees of anisotropy. We choose this example because it models a weakly bound, strongly anisotropic molecule that exhibits large amplitude motion. Thus, there is no obvious “best” zero-order approximation to describe its states. The vibrational-rotational motion was solved variationally by using a body-fixed Hamiltonian. Rather than using harmonic or Morse oscillator basis functions for the van der Waals stretching coordinate, we used basis functions constructed from products of polynomials and the ground-state wave function of a LennardJones potential. The total wave function was then transformed and expressed in terms of natural orbitals. As expected, the natural orbital expansion of the wave function converged much more quickly than the original expansion, thus allowing for easier interpretation of the wave function. In order to better understand the results of the natural orbital analysis we also made use of the harmonic limit of the wave equation by relating the nodal behavior in this limit to that of the natural basis function expansion.
(1) D. H. Levy, Annu. Reu. Phys. Chem., 31, 197 (1980). (2) G. D. Carney, L. L. Sprandel, and C. W. Kern, Adu. Chem. Phys.. 37, 305 (1978). (3) R. J. LeRoy and J. S . Carley, Ado. Chem. Phys., 42, 353 (1980); see also the original work: C. F. Curtiss, J. 0. Hirschfelder, and F. T. Adler, J . Chem. Phys., 18, 1638 (1950). (4) A. M. Dunker and R. G. Gordon, J . Chem. Phys., 64, 354 (1976). (5) J. M. Hutson and B. J. Howard, Mol. Phys., 41, 1123 (1980). (6) R. J. LeRoy and J. Van Kranendonk, J . Chem. Phys., 61,4750 (1974). (7) P. 0. Lowdin, Phys. Reo., 97, 1474 (1955). (8) E. R. Davidson, Reo. Mod. Phys., 44, 451 (1972). (9) E. R. Davidson, “Reduced Density Matrices in Quantum Chemistry”, Academic Press, New York, 1976. (IO) A. E. Reed and F. Weinhold, J . Chem. Phys., 78, 4066 (1983). (1 1) C. Edmiston and K. Ruedenberg, Reo. Mod. Phys., 35, 457 (1963). (12) P. D. Robinson, J . Chem. Phys., 66, 3307 (1977). (13) D. M. Bishop and L. M. Cheung, Phys. Rec. A, 16, 640 (1977). (14) D. M. Bishop and L. M. Cheung, In?.J . Quanrum Chem., 15, 517 (1979). (15) R. M. Stratt, N. C. Handy, and W. H. Miller, J . Chem. Phys.. 71, 3311 (1979). (16) M. L. Sage and J. A. Williams 111, J . Chem. Phys.. 78. 1348 (1983).
0 1985 American Chemical Society
The Journal of Physical Chemistry, Vol. 89, No. 21, 1985 4465
Vibration-Rotation Wave Functions Before specializing the discussion to the analysis of wave functions, in the next section we review the general formulation of the wave equation for triatomic van der Waals molecules in body-fixed coordinates. In section 3 we review natural orbital analysis and present a specific application of it to obtain natural orbitals of the vibrational degree of freedom. A model potential based on the molecular parameters of He-LiH is introduced in section 4 along with details of the calculations, and in section 5 we present the results of our calculations and the analysis of the wave functions. A brief summary concludes the paper.
the I;*are the spherical harmonics, and the AiujSujn are the expansion coefficients to be determined. Substituting the expansion 7 into the set of coupled differential eq 2 reduces to a secular equation
where
2. The Wave Function Our formulation of the wave equation is essentially identical with that of the body-fixed Hamiltonian customarily used for triatomic van der Waals molecule^,^ but for completeness we will review the important equations. We have used a standard coordicate system in which 7 is the vector from nucleus A to nucleus B, R is the vector from the center of mass of the diatom, AB, to nucleus C , and 0 is the angle between R and 7. In this coordinate system stationary-state eigenfunctions of the total angular momentum have the form
where J is the total angular momentum, M is its space-fixed z component, and d,, is a representation coefficient of the rotation from space-fixed to body-fixed coordinates specified by the Euler angles a,6, and y. The internal wave functions xi satisfy the following set of coupled differential equations3J7
h 2X+( J,Q)j-(n) Xi+l
hZX-(J,Q)j+(n) 2pR2
-
h2[J(J 2pR2
+ 1) + j ( j + 1) - 2Q2] + + Vd(r)
h2j(i 1) 2pdr2
d2 + 2kd dr2
9JMp
(2)
In the derivation of these equatkns it has been assumed that the
g,,(r) d r
+
For a fixed basis set, eq 8 and 9 are the working equations which determine the wave function. There is, however, one additional symmetry that can be recognized from the outset, and that is the invariance of the total Hamiltonian to inversion of the coordinate axes. This symmetry allows the internal wave function to be grouped into parity eigenstates, and the total wave function with parityp 07 = 0 is even, p = 1 is odd) is J
xi-1 = 0
1
= C C*[XJnDJn, R=p
+ (-1)P
X!&*M]
(10)
where
c*=1/2, n = o
body-fixed z axis *coincides with R, and that r‘ lies in the body-fixed x z plane. The Hd is the diatomic Hamiltonian
= 1/2’/2, Q
>0
(11)
If eq 10 and 11 are used, the secular equation may be rewritten as where
pd
is the reduced mass of the diatom
J
x x (JMpnvjnlHIJMpn‘v3‘n‘)Ai~u~n, = EJpA${,R
R’=p
is the angular momentum of the diatom expressed in kody-fixed coordinates and tl is the quantum number of J, and j , . In the coupling terms X*(J,Q) = [J(J
+ 1) - n(n i 1 ) ] 1 / 2
(5)
where (JMpnujl?lHlJMpn’u[iQ’) = 0
for p = 1 and Q or
where thef, and g, are appropriate orthonormal basis functions, (17) M. Shapiro and G. G. Baht-Kurti, J . Chem. Phys., 71, 1461 (1979).
n’ = 0 (13a)
= 21i2(JMnvjnlHIJMn’u’Q’) for p = 0 and (Q,Q’) = (0,l) or (1,O) (1%) = (JMnujQlHIJMn’u’j’Q’) otherwise
Solutions to eq 2-6 are obtained by either the method of closecoupled equations”6 or the variational method.2 For the calculations presented in section 5 , we have used the variational method. In the variational approach, the internal wave functions x i are represented in terms of a basis set2
(12)
fl’b$#
(13c)
and the Ai{,* are the expansion coefficients which are determined in our calculations. These coefficients define the total wave function, and therefore, as we show in the next section, the Ai{,* play an important role in the natural orbital analysis.
3. Natural Orbital Analysis In electronic structure theory the first-order reduced density matrix is defined as7-9 y( 1,l’) = N s 9 * ( 123 ...N ) 9(1’23...N)d2d3 ...d N (14)
where the numerals refer to electron coordinates and the natural orbitals are defined as the eigenfunctions of
4466
The Journal of Physical Chemistry, Vol. 89, No. 21, 1985 S y ( 1 , l ’ ) q(1’)dl’ = Xiui(l)
(15)
For nuclear motion we make use of the density matrix for a specific coordinate, and the eigenfunctions of the density matrix are the natural orbitals for the coordinate. A better terminology in the general case is perhaps “natural basis functions” since “orbitals” tend to imply electronic motion. AnotherI6 suggestion has been “natural local modes”. One must keep in mind that the results of a natural orbital analysis, when applied to vibrational degrees of freedom, are highly dependent on the particular coordinates used. For example, if we used a space-fixed axis system and a different set of internal coordinates, there would be no simple relation between the natural basis functions obtained and those which we discuss below. The analysis is most useful when the coordinates are those in which the Hamiltonian is nearly separable. Thus, a local mode analysis may not be the best for a particular molecule. With this comment, we shall revert to the use of natural orbitals terminology. The derivation of natural orbitals for the van der Waals bond length R begins with the density matrix p(R,R3 = l*JP*(R,r,8,a,B,~) X *JP(R’,r,8,a,&y) d a sin B dB d y sin 8 d8 r2 d r = J
c Pn(R,R?
R =p
(16)
where we have decomposed the density matrix into its components for each Q pn(R,R? = Sx#*(R,r,8) x$’(R’,r,O) sin 8 d0 r2 d r (17) This decomposition is always valid, and for strongly anisotropic potentials, when il is nearly a good quantum number, only one pn is dominant. Thus, we define natural orbitals for each Q,using PR. When the x f ’ s are expressed in terms of a basis set, as in eq 7, then Po(R,R? = eB;fl,f”*(R)
fn@?
(18)
nn‘
where
Bin, = ZAi$;R Ai{,R
(19)
LJ
The eigenvalue equation defining the natural basis functions for the R coordinate is BOQR = QOAn
(20)
where the elements of Bn are the Bf,,,, and the natural basis functions are Ffla(R) = q A R ) n
QL
(21)
The eigenvalues, AR, are the “occupation numbers” or relative importance of each natural orbital, and normalization of the wave function sets C T r An = 1 R
(22)
Likewise, the density matrix of 0 for a given Q is pn(8,8’) =
1 x#*(R,r,B) x$’(R,r,8’)RZ dR ? d r = CGt ~ J n ( 8 , O Y,,n(e’,o) )
and GuR(r)and SiTuare the natural functions and transformation coefficients for the r coordinate, which we are not treating explicitly. We will use the form of the internal wave function given by eq 27 in our analysis of the wave function. Although there are just as many terms in this transformed expansion as in the original, as we will see, only a few are dominant, so a significant simplification in the wave function results. 4. Calculations
In this section we present a model potential for anisotropic van der Waals molecules which allows us to examine bound states for both the linear and T-shaped equilibrium configurations along with various degrees of anisotropy. The masses and other molecular parameters used are based on those of He-LiH, the simplest polar van der Waals molecule. This van der Waals system, in principle, provides an example for which both theory and experiment can be accurately applied to determine bound states, although experimental results are currently limited to high-energy scattering A. Model Potential. Extensive ab initio calculations were performed on the He-LiH system by S i 1 ~ e r . IEight ~ Slater-type basis functions were used, and many-body perturbation theory was applied through third order in the correlation energy. The ab initio evaluation of the surface was performed at 144 points required by a numerical quadrature scheme which yielded an analytic fit to the potential energy surface. This surface has a linear equilibrium geometry with the He located near Li at a distance of 4.56 bohr from the center of mass of LiH (see Figure la). The well depth is 15.3 meV or about 123 cm-I, and no potential well was found at the H end (see Figure la). Rather than using the analytic form of the potential surface of He-LiH given by Silver, we adopted a highly idealized potential which can be easily modified to study different degrees of anisotropy. We felt free to modify our potential, since our main purpose was not to calculate highly accurate energy levels, but to demonstrate the utility of natural orbital analysis for giving insight into the vibrational motion of a rotating molecule. We averaged over the vibrational motion of the diatom, thus leaving two internal degrees of freedom, R and 8. The interaction potential U(R,8) was represented by a modified Lennard-Jones 6-12 potential
(
;)12
- (1
+ X cos 8 - A)
for the linear equilibrium configuration, and
c”,, = y
;
R
Ai$,%
(24)
U(R,8) = 4c[
Then the natural orbitals for the 8 coordinate are
where
where
U(R,8) = 4c[
where
=
P P = Pro (26) and the occupation numbers are given by I’*. The overbar is used in the definition of the natural functions for the 8 coordinate to denote that these functions are averages of the spherical harmonics. The label j in eq 25 is not an angular momentum quantum number, however. The natural orbitals for the r coordinate are defined in a similar fashion; however, since we will use only a single basis function for the r coordinate in our work, we will not give the equations here. With the transformations given by eq 21 and 25, the internal wave functions can now be written in terms of the natural orbitals
(23)
JJ
qn(‘)
Holmer and Certain
C?‘n(8,0) FJ J,
(25)
(
;)I2
- (1
+ X sin 8 - A)
($1 ($1 -
(29)
-
(30)
(18) P. J. Dagdigian and B. E. Wilcomb, J . Chem. Phys., 72,6462 (1980). (19) D. M. Silver, J . Chem. Phys., 72, 6445 (1980). (20) E. F. Jendrek and M. H. Alexander, J . Chem. Phys., 72,6452 (1980).
The Journal of Physical Chemistry, Vol. 89, No. 21, 1985 4461
Vibration-Rotation Wave Functions
mH'
6
mH
3 0
-3
mH -6
i
I
4
-3
0
3
mH
s
0
-6
-3
0
3
6
Figure 1. Comparison of potential energy contour diagrams of the He-LiH interaction surface: (a) analytic representation of Silver,I9 and (b) eq 29 with X = 1.9. Successive contours differ by a factor of lo'/* in energy (expressed in millihartree), and distances shown are in bohr units.
Rather than usingf,'s obtained by Gram-Schmidt orthogonalizing the fn's, we used the eigenvectors of Ho obtained from
TABLE I: Physical Constant Used in the Model Calculations m (4He) 4.00260 amu m ('LiH) 8.023825 amu Re = 2'/'0 2.413 8, Bo ('LiH) 7.4065 cm-' c123.4 cm-'
Ho W = S WE0 where
for the T-shaped configuration, and the value of A governs the degree of anisotropy. Values of the physical constants used in our calculations are given in Table I. B. Basis Set. In all of our calculations we made the approximation of using only the ground-state wave function go(r) for the LiH vibrational motion. Thus, in effect, we limited ourselves to a rigid diatomic rotor. Averaging over the diatomic vibrational motion, we obtained
and
where eo and Bo are the LiH's ground-state vibrational energy and rotational constant, respectively. Rather than using harmonic or Morse oscillator functions for the f n ( R ) ,we determined them by the method given by Sitz and Yaris2I In this method the basis functions are obtained by orthogonalizing a set of trial functions (Pn(R) do(R)] where the Pn(R) are functions of R, and do(R) is the ground-state wave function of some reference Hamiltonian, Ho(R). Matrix elements for the R coordinate may be calculated from integrals of the type I d o ( R ) h(R) do(R) dR, where h(R) is some function of R . In our application we have chosen Ho(R) = --h 2 d2 211 dR2
+ 4€[ (;)I2
-
($1
(33)
since this incorporates the R-dependent part of the wave equation 2 which is independent of 0, A, J , and Q . The ground state, do(R), was determined by the Cooley methodz2on a 501 point grid from 0.50 to 5 . 0 ~ .We chose R" for the Pn(R), and therefore the set of product functions was
(21) P.Sitz and R.Yaris, J. Chem. Phys., 49, 3546 (1968). (22) J. W. Cooley, Math. Comp., 15, 363 (1961).
(35)
and
and then fn(R) = Cfn,(R) n'
Wdn
(38)
This orthogonal set of fn's diagonalizes Ho and allows matrix elements of this part of the total Hamiltonian to be replaced with Eom6nd. The above procedure was done only once for all of the calculations, and it simplified the matrix elements of the total Hamiltonian. There were, however, problems with linear dependence of thefis, due to our choice of Rn for the P i s , which limited the number of f i s to seven. The total Hamiltonian matrix was constructed by using seven f,(R), and spherical harmonics, qn(O,O), with j ranging from Q to 8, and Q's from p to J. All possible products of the fn's and 5s ; were taken and resulted in 63- to 168-term expansions for the examples studied here. 5. Results and Discussion In this section we apply the methods presented in the previous sections to three specific examples: two linear equilibrium geometries, but with different degrees of anisotropy, one of which mimics the potential surface for He-LiH given by Silver,I9 and one with a T-shaped equilibrium geometry. We also present a harmonic analysis of the wave equation in order to interpret the natural orbital results. A. Linear Configuration with Strong Anisotropy. With a value of 1.9 for the anisotropy parameter A, the model potential, eq 29, compares reasonably with the representation of the He-LiH potential surface given by SilverI9(see Figure l), and it also retains three bound vibration-rotation states whose energies are given in Table 11. The energies listed in Table I1 are accurate to about 0.01 cm-'. These eigenstates are bound by only a few percent of the total well depth and are very sensitive to the value of A.
4468
The Journal of Physical Chemistry, Vol. 89, No. 21, 1985
TABLE 11: Energy Levels for the Model Potential Given by Eq 29 with X = 1.9 J 0 1
parity even even
2
even
energy, cm-'
Holmer and Certain TABLE 111: Ten Largest Coefficients of the Basis Set Expansion (Eq 7) for the Model Potential Given by Eq 29 with J = 0, p = 0, and h = 1.9O ~~
-5.375 -4.006 -1.302
n
j
0
1 0 2 0 1 3 2 4 1 1
1
A$* 0.592 86 0.478 08 0.413 10 0.32001 0.271 71 0.210 00 0.095 02 0.084 99 0.060 98 -0.049 25
OThe total number of coefficients is 63.
62
U
? 0
TABLE I V Five Largest Coefficients of the Natural Orbital Expansion (Eq 27) for the Model Potential Given by Eq 29 with J = 0, p = 0, and h = 1.9'
G
-1
n
j
a
.A:?n
0 1 2
0 1 2 3 4
0 0 0 0 0
-0.987 05 -0.1 55 95 -0.036 50 0.008 93 0.001 79
3 4
1
2
'The total number of coefficients is 63.
3
R
Figure 2. First four natural orbitals of the R coordinate for the eigenstate with J = 0 and p = 0 calculated by using the model potential given in eq 29 with X = 1.9.
TABLE V Convergence of the Truncated Natural Orbital Expansion (Eq 27) for the Model Potential Given by Eq 29 with J = 0, p = 0, and X = 1.9" no. of terms E , cm-' 63
no. of terms
6
1 2 3 4 5
1-
U
0
I
R 0 0 0 0 0 0 0 0 0 0
s
0-
-5.375
Complete Wave Function (R) ( u ) (R2) (a2) (cos e ) (P2(c0se)) -1.3973
-2.0045
0.8254
0.5763
Truncated Natural Orbital Expansion AE A(R) A(R2) A(cos 8 ) A(P2(cose ) ) 3.331 0.343 0.032 0.002 0.000
-0.0058 -0.005 0.0000 0.0000 0.0000
-0.0213 -0.0022
-0.0002 0.0000 0.0000
0.0107 0.0011 0.0001 0.0000 0.0000
0.0142 0.0006 0.0000 0.0000 0.0000
"The quantities given are the errors (A) in each of the expectation values.
2
1
3
9 Figure 3. First four natural orbitals of the 0 coordinate for the eigenstate with J = 0 and p = 0 calculated by using the model potential given in eq 29 with X = 1.9.
Thus, there is a distinct possibility that the real He-LiH system has no bound vibration-rotation states, but this requires careful investigation with a more accurate potential and a larger basis set. We calculated the natural basis functions, as outlined on section 3, for the ground-state J = 0 wave function. The first four natural functions of R and 0, ordered by their occupation numbers, are shown in Figures 2 and 3. We find that the natural orbitals of R are reminiscent of the eigenfunctions of a Lennard-Jones potential, yet the excited states of the latter emphasize the outer turning point region much more than the corresponding natural orbitals. The natural orbitals of 0 make a transition from harmonic oscillator-like wave functions (j= 0) to hindered rotor wave functions (j= 3). The internal wave function, expressed in terms of the natural orbital expansion, is summarized in Table IV, and can be com-
pared with the ten largest terms of the original basis set expansion (eq 7) given in Table 111. Many of the coefficients of the original basis set expansion are the same size and the magnitudes of the coefficients do not drop off very quickly. On the other hand, there is one dominant term in the natural orbital expansion and the magnitudes of the coefficients of the subsequent terms decrease rapidly. The natural expansion also shows no "cross terms", that is, terms such as (n = 0 , j = l), (n = 2 , j = l), etc. This is a special property of the natural orbital expansion of two coordinate wave functions (cf. two-electron wave functions in electronic structure theory9). As a measure of the rate of convergence of the natural expansion, the expansion was truncated, renormalized, and used to calculate molecular properties. This is shown in Table V, where, for example, truncation to three terms means that the first three given by Table IV were used. We find that, with each additional term included, the error in all the expectation values is reduced by an order of magnitude. Interestingly, the relative error in the binding energy is quite large for the one term wave function, whereas the other expectation values are corrected to within 2%. Contour plots of the individual terms of the natural orbital expansion, i.e. (39)
are instructive, and these are given in Figures 4-6. Positive contour values are represented by solid curves, negative by dashed curves, and zero by long dashed lines. The dominant term of the
The Journal of Physical Chemistry, Vol. 89, No. 21, 1985 4469
Vibration-Rotation Wave Functions
2.5
I-
I
n
D
v
I
I
e 1.5
1.0
I 45
90
c
I
90
45
135
e
e (degrees) Figure 4. Contours of the individual term of the natural orbital expansion (eq 39) with n = 0, j = 0, Q = 0, and p = 0 calculated by the model potential given by eq 29 with X = 1.9. Contour values are equally spaced 0.5.
I
135
[degrees)
Figure 6. Contours of the individual term of the natural orbital expansion (eq 39) with n = 2, j = 2, Q = 0, J = 0, and p = 0 calculated by using the model potential given by eq 29 with X = 1.9. Contour values are equally spaced by 0.02. TABLE VI: Seven Largest Coefficients of the Natural Orbital Expansion (Eq 27) for the Model Potential given by Eq 29 with J = 2. D = 0. and X = 0.4” 0 0 1 1 0 2 2
n
0
v
P
1 0 2 1 2 3 2
1 0 1 0
2 1 0
0.920 80 -0.358 78 -0.105 62 -0.10067 0.044 08 0.010 54 -0.006 59
“The total number of coefficients is 168. TABLE VII: Convergence of the Truncated Natural Orbital Expansion (Eq 27) for the Model Potential Given by Eq 29 with J = 2, p = 0,and X = 0.4” 45
90
135
8 [degrees)
Figure 5. Contours of the individual term of the natural orbital expansion (eq 39) with n = 1, j = 1, Q = 0, J = 0, and p = 0 calculated by using the model potential given by eq 29 with X = 1.9. Contour values are equally spaced by 0.1.
expansion (n = 0, j = 0, R = 0) is simply a product of ground-state oscillator functions. Each additional term has one more node in both the R and 8 directions than the previous term. One can superimpose the graphs to determine how each correction term adds to or subtracts from the sum of all the previous terms (note, however, that the contour value spacings are different). The general trend in this case is to add to the wave function along a line from small R and small 8 to large R and large 8. This is due to the shape of the potential and is discussed below. The two other eigenstates listed in Table I1 are the lowest vibrational states with J = 1 and 2. Natural analysis of their wave functions gives results similar to those for the J = 0 ground state given above. There are, however, terms in the natural expansion with R # 0. Rather than describing these states farther, we move on to a less anisotropic potential and study an example of excited vibrational motion with J = 2. B . Linear Configuration with Medium Anisotropy. To illustrate the results of natural orbital analysis when the distinction between different R components of the wave function is important, we decreased the value of X to 0.4 and considered an excited bending state with J = 2 and p = 0. The most important terms of the natural orbital expansion are summarized in Table VI. The dominant term in the expansion is an R = 1 component, which
no. of terms E , cm-I 168
no. of terms 1 2 3 4 5 6 7
-5.890
Complete Wave Function ( R ) (a) (R2) ( u 2 ) (cos 0 ) 1.3778
1.9457
0.3463
(P2(cose ) ) -0.0567
Truncated Natural Orbital Expansion AE A(R) A ( R 2 ) A(cos 0 ) A(P,(cos 8 ) ) 2.096 1.386 0.407 0.103 0.034 0.012 0.003
-0.0293 -0.0013 -0.0005 0.0001 0.0000 0.0000 0.0000
-0.0936 -0.0057 -0.0027 0.0003 -0.0001 -0.0001 0.0000
0.0674 -0.0015 -0.0046 0.0002 0.0001 0.0000 0.0000
0.OE -0.0084 -0.0073 0.0002 -0.0001 0.0000 0.0000
”The quantities given are the errors (A) in each of the expectation values. indicates that the bending vibration has one quantum of vibrational angular momentum about the body-fixed z axis. We find that the first correction term to the natural expansion is an fl = 0 component, and there is also a minor contribution from a term with R = 2. The convergence of the natural expansion is presented in Table VII, which shows that the errors in the various expectation values decrease rapidly but nonuniformly. The error in (cos 0 ) actually increases when going from a two- to a three-term wave function. Contour plots of the individual terms of the natural expansion (eq 39) are shown in Figures 7-10. As mentioned above, the most significant term is an excited bending mode with one quantum of vibrational angular momentum. The second term is also an excited bending state, but with Q = 0. Both terms correspond to ground-state oscillators along the van der Waals
4470 The Journal of Physical Chemistry, Vol. 89, No. 21, 1985
-
Holmer and Certain
2.0-
0 v
c: 1.5'
1.0-
I 45
90
135
45
90
8 (degrees)
135
8 (degrees)
Figure 7. Contours of the individual term of the natural orbital expansion (eq 39) with n = 0, j = 1, Q = 1, J = 2, and p = 0 calculated by using the model potential given by eq 29 with X = 0.4. Contour values are equally spaced by 0.2.
Fngure 9. Contours of the individual term of the natural orbital expansion (eq 39) with n = 1, j = 2, Q = 1, J = 2, and p = 0 calculated by using the model potential given by eq 29 with X = 0.4. Contour values are
equally spaced by 0.05. J
1
I
(e-
2.0'
I I 1 I -----+--_---I
h
.
0
Y
\
1.5'
I I 1 I I
t
i
2.51
/I c
1.0-
1 ' .
t
I I
45
I
1
90
.-----_ *-----
I I I
135
90
45
8 [degrees)
4
I 135
8 (degrees)
Figure 8. Contours of the individual term of the natural orbital expansion (eq 39) with n = 0, j = 0, 0 = 0, J = 2, and p = 0 calculated by using
Figure 10. Contours of the individual term of the natural orbital expansion (eq 39) with n = 1, j = 1, Q = 0, J = 2, and p = 0 calculated
the model potential given by eq 29 with X = 0.4. Contour values are equally spaced by 0.1.
by using the model potential given by eq 29 with X = 0.4. Contour values are equally spaced by 0.05.
TABLE VIII: Four Largest Coefficients of the Natural Orbital Expansion (Eq 27) for the Model Potential Given by Eq 30 with J = 1. D = 0. and X = 1.0" n j R .xfn
TABLE I X Convergence of the Truncated Natural Orbital Expansion (Eq 27) for the Model Potential Given by Eq 30 with J = 1, p = 0, and X = 1.W
0 1 0 2
0 1 1 2
0 0 1 0
-0.997 48 -0.066 52 0.024 31 -0.004 60
"The total number of coefficients is 119. stretching coordinate, R. The third term in the natural expansion (Figure 9) is again an Q = 1 component which corrects the most significant term (Figure 7) in a manner analogous to the first example, adding to the wave function along a line from small R and small 6. The fourth natural term (Figure 10) is an Q = 0 term, just as the second term. We see that these two terms are of comparable magnitude, yet emphasize very different regions in the R,B-space, which is indicative of nonseparability (at least in this coordinate system). C. T-Shaped Configuration with Strong Anisotropy. To illustrate the results of natural orbital analysis for a T-shaped molecule, we used the model potential, eq 30, with X = 1.0. The
Complete Wave Function ( R ) (a) (R2)(a2)
no. of terms
E, cm-'
119
-43.446
no. of terms 1 2 3 4
1.2755
1.6496
(P,(cos 8 ) ) -0.2590
Truncated Natural Orbital Expansion AE A(R) A ( R 2 ) A(P2(cos 8 ) ) 0.747 0.047 0.009 0.000
-0.0003 0.0001 0.0000 0.0000
-0.001 1 0.0001 0.0000 0.0000
-0.0028 -0.0002 0.0000 0.0000
"The quantities given are the errors (A) in each of the expectation values. four largest terms in the natural orbital expansion for the lowest vibrational state with J = 1 and p = 0 are given in Table VIII. Here we find a slight contribution from an D = 1 component, but the amount is much less than the contributions of Q = 0 components in the previous example. Convergence of various expectation values is illustrated in Table IX. The errors drop by
The Journal of Physical Chemistry, Vol. 89, No. 21, 1985 4471
Vibration-Rotation Wave Functions
I I I I I I
I I I I I I I I I I I
h
0
v
e
45
90
135
45
8 (degrees)
2*sl
4s
I I
I
I
I I
I
I 90
Figure 13. Contours of the individual term of the natural orbital expansion (eq 39) with n = 0, j = 1, s2 = 1, J = 1, and p = 0 calculated by using the model potential given by eq 30 with X = 1.0. Contour values are equally spaced by 0.01.
the R coordinate, but a change of one along the 0 coordinate. In order to gain insight into this behavior, we examine the smallamplitude limit of the wave equation and relate the various terms of the wave equation to harmonic oscillator operators. Even though our model systems exhibit wide amplitude motion, we find that the symmetry and nodal properties of the harmonic limit explain the behavior of the natural orbitals expansion. The small amplitude limit of the wave equation 2 about R = Re and 8 = 0 (for the linear configuration) is
I I
I
135
8 (degrees)
Figure 11. Contours of the individual term of the natural orbital expansion (eq 39) with n = 0, j = 0, s2 = 0, J = 1, and p = 0 calculated by using the model potential given by eq 30 with X = 1.0. Contour values are equally spaced by 0.5. I
90
c
a ~ + i -ae - - ) oR X ~ + ~ -
135
8 (degrees)
Figure 12. Contours of the individual term of the natural orbital expansion (eq 39) with n = 1 , j = 1, s2 = 0, J = 1, and p = 0 calculated by using the model potential given by eq 30 with X = 1.0. Contour values are equally spaced by 0.025.
about an order of magnitude with each additional term in the expansion, yet again we find that the relative error in the binding energy is larger than the relative error in the other expectation values. The average of (cos 0 ) is not listed because it is zero. Contour plots of the individual terms of the natural expansion (eq 39) are given in Figures 11-13. The most significant term is a product of ground-state oscillator functions centered about 0 = 90'. The second term corrects the first by adding to regions 0 ' 1 to large R and along the lines from small R and small 10 - 9 large 10 - 90'1. The third term is an 52 = 1 component with no nodes in the R direction but one in the 0 direction. D. Harmonic Analysis. As we have seen, natural orbital analysis allows the wave function to be. expressed as a dominant term plus several correction terms. However, natural orbital analysis does not give an explanation for the magnitude and ordering of the correction terms. For the sequence of natural orbitals with a given value of Q , we expect that each new term in the expansion would have one (or two if dictated by symmetry) more nodes along each of the two coordinates. This is exactly what we have seen in the previous examples. However, the major correction term with a value of 52 different than that of the most significant term shows no change in the number of nodes along
The wave equation is thus seen to reduce to a simple harmonic oscillator in the R direction and the radial component of a twofold degenerate oscillator in the 0 direction, plus potential energy and angular momentum coupling terms. The oscillator eigenfunctions with the same 52 are coupled by a cubic term from the interaction potential, whereas the kinetic terms couple oscillator wave functions with different values of 52. These kinetic terms are just linear combinations of raising and lowering operators2j
a ao o - - a- - = -1;2 + 1 - - - =52- - 1
ao
o
it
- ill+ 2
iv - it' 2
(41)
where t@(V,Q)a @(V+l,Q+l) tt@(V,52) = @(V-l$-l) v@(V,Q) a @(V+l,Q-l) ?$@(V,Q)
a
@(V-l,Q+l)
(42)
(23) J. D. Louck and W. H. Shaffer, J . Mol. Spectrosc., 4, 285 (1960).
4472
The Journal of Physical Chemistry, Vol. 89, No. 21, 1985
and where V and R are the energy (or radial) and angular momentum quantum numbers of the twofold degenerate harmonic oscillator eigenfunctions, @( V,R). Thus, in first-order perturbation theory, two types of coupling terms between the zeroth order harmonic oscillator eigenfunctions exist: potential coupling, which couples harmonic oscillator functions that have the same R and differ by one in the R quantum number (or number of nodes) and zero or two in the 6 quantum number; and kinetic coupling, which couples states which differ by one in R and one in the 6 quantum number with no change in the R quantum number. For the T-shaped geometry, the small amplitude limit of the wave equation is not a twofold degenerate oscillator in the 8 direction; however, the kinetic coupling has the same effect on the zeroth-order harmonic oscillator eigenfunctions as the kinetic coupling for the linear geometry. Returning to the results of the natural orbital analysis, the nodal structures of the most important terms in the natural expansion (for a given R) differ in the same manner as the nodal structures of the harmonic zero-order states coupled by the potential coupling term. This suggests that these terms are corrections for nonseparability in the potential. In all three examples we found that these correction terms added to the most significant term along a line from small R and small 16 - Oesl to large R and large (66,1 which is required by the direction of the skew in the potential. For example, for linear geometry the minimum of U(R,B)for a given 6 is
This shows that, as 6 deviates from its equilibrium position, R&(6) becomes larger. Thus, the correction terms tend to add to the wave function along this line. The first correction term in the natural expansion with a value of R different than that of the most significant term has a nodal structure which differs from that of the most significant term in the same way as harmonic zero-order states coupled by the kinetic coupling terms. Thus, these correction terms (e.g., the (n = 0, j = 0, R = 0) term in Table VI) arise from nonseparability introduced by the kinetic coupling or coupling of overall rotation with the internal motion. These kinetic coupling terms in the wave eauation are ODerators of both R and 6, but R appears only as l j R 2 . Near R ’= Re, 1 / R 2is a slowly varying function, suggesting the reason why there is no nodal change in the R direction for these types of natural expansion correction terms.
Holmer and Certain
6. Summary We have presented an application of natural orbital analysis to vibration-rotation wave functions. The expansion of the total wave function when expressed in terms of the natural orbitals converged quickly and facilitated interpretation of the various components of the wave function. For the specific model potentials used in our calculations we related the contributions of potential and kinetic coupling to the nodal structure and spatial distributions of the correction terms. In the first example all corrections were due to potential coupling, since there was no total angular momentum (J = 0). The second and third examples illustrated the interplay between potential and kinetic coupling: for the second example (J = 2, medium potential anisotropy) the largest correction was due to kinetic coupling, whereas for the third example (J = 1, strong potential anisotropy), potential coupling gave rise to the largest correction term. For the model potential used in our work, R was nearly a good quantum number, and therefore it was useful to calculate natural basis functions for each R since. only one po was dominant. Weakly anisotropic systems are more appropriately described by spacefixed coordinates, and it may be more useful for these systems to base the analysis on natural functions of (O1,4Jand (62,42), where Bi and $J~are the polar angles of R and 7. Another possibility is to find natural functions defined by diagonalizing p(B,a,P,y;O’,a’,pI,y’) or p(61,41,62,42;61’,41’,6~,4~). These functions make y/,(BZ,+J a smooth transition from the spaced-fixed Ylm,(O1,+,) to the body-fixed qn(6,0) d,,(a,P,y), but are functions of four variables and cannot be visualized easily. As mentioned earlier, the results of a natural orbital analysis, when applied to vibrational degrees of freedom, are highly dependent on the particular coordinates used. In practice, one should use coordinates which make the Hamiltonian almost separable. One possibility is to optimize the coordinates self-consistently as suggested by Thompson and T r ~ h l a r . Another ~~ possibility is to choose coordinates which maximize the contribution of the dominant term in the natural expansion. These optimized coordinates could be used to study the transition from normal mode to local mode vibration^.^^ Acknowledgment. This research was supported by NSF Grant NO. CHE 8 1-07628. (24) T.C. Thompson and D. G. Truhlar, J . Chem. Phys., 77,3031(1982). (25) R.T.Lawton and M. S.Child. Mol. Phvs.. 40.773 119801: 44. 709 (i&ij.
