Empirical potential energy surface for Ar·SH/D and Kr·SH/D

Empirical potential energy surface for Ar·SH/D and Kr·SH/D Prakashan P. Korambath, Xudong T. Wu, Edward F. Hayes, Christopher C. Carter and Terry A. Miller Department of Chemistry The Ohio State University Columbus Ohio 43210

Abstract Experimental data from vibrationally and rotationally resolved laser induced fluorescence experiments have been used to produce potential energy sure 2 Σ+ states of the Ar·SH and Kr·SH van der faces (PES) for the excited A Waals complexes. This was done using a potential energy functional form first suggested by Bowman and co-workers1 for Ar·OH/D. A discrete variable representation (DVR) of the vibration–rotation Hamiltonian was used in combination with the Implicitly Restarted Lanczos Method and Sequential Diagonalization Truncation (SDT) of the DVR Hamiltonian. This approach takes advantage of the sparseness of the DVR Hamiltonian and the reduced order of the SDT representation. This combination of methods greatly reduces the amount of computational time needed to determine the eigenvalues of interest. This is important for the determination of the PES that results from minimizing the difference between the experimental and theoretically predicted values for the vibronic energy levels and their corresponding rotational constants. In addition this procedure was helpful in assigning the absolute vibrational quantum numbers for the deuterated species for which less experimental data was available. Plots of the calculated wavefunctions cor-

1

responding to various experimental vibronic bands indicate that these states sample regions of the PES from 0 degrees, where the hydrogen atom is closest to the rare gas atom, to approximately the saddle point, near the T-shaped configuration. As a result this region of the surface is determined accurately whereas the region of the PES around 180 degrees, corresponding to the sulfur atom being closest to the rare gas atom, is determined only qualitatively.

2

I. INTRODUCTION

We have recently undertaken a comprehensive experimental and theoretical effort to characterize the R·SH/D (R=Ne, Ar, and Kr) complexes. Experimentally we have assigned, via heavy atom isotope shifts and other techniques, a large number of vibronic bands of e 2 Π3/2 (0,00 ,0) for these complexes.2 Higher e 2 Σ+ (0,vk ,vs ) – X the electronic transition A b resolution spectra3 of these bands have complemented the vibrational assignments. In addition, they have yielded highly precise values of the rotational, fine, and hyperfine constants of these levels. A complementary resolved emission study,4 albeit of lower resolution, has e state of the complexes. yielded considerable information about the X Our recent experimental work on the R·SH species in turn complements the existing, significant database on the R·OH/D (R=Ne, Ar, and Kr) complexes.5 It is probably safe to say that taken together the R·XH/D (X=O, S) are the best experimentally characterized complexes known. The R·XH/D complexes thus represent an extremely fertile area for theoretical approaches to produce reliable potential energy surfaces (PES). Indeed the data available for these complexes hold promise for the most reliable PES, based upon experiment, and also a stringent test of theoretical approaches. Several efforts along these lines have already been published5–11 for the lighter R·OH complexes. However, similar work on the thio analogues has thus far been lacking. The specific purpose of this paper is to derive accurate PES for the moderately strongly e 2 Σ+ states of Ar·SH and Kr·SH based on the recently reported experimental bonded A results. The functional form of the potential energy surface is taken from the earlier work of Bowman and co-workers.1 This form, which was first used for Ar·OH, has several adjustable parameters that govern the potential energy surface. For Ar·SH, and Kr·SH, we determined these parameters by minimizing the root–mean–square difference between the calculated and experimental results. To calculate the rovibronic energy levels efficiently, we use a discrete variable representation (DVR) approach originally introduced by Light and co3

workers.12 The computational aspects of this original approach have been enhanced in three major ways. First, we use the Implicitly Restarted Lanczos Method (IRLM) of Sorensen et al.13 to determine the eigenpairs of interest. Second, we use a Sequential Diagonalization Truncation (SDT) approach introduced by Light and coworkers 14 to reduce the order of the matrix representation. Third, we apply a Chebychev polynomial preconditioning to speed up the convergence of the IRLM.15,11 Friesner et al.16 and Bramley and Carrington17 have both independently reported on the value of the implicit SDT transformation that was introduced by Pendergast et al.15

II. THEORETICAL APPROACH

A. Potential Energy Surface

The global potential for the triatomic system A–BC is based on a set of three one dimensional radial potentials that are combined using switching functions. The two minima, at 0 degrees and at 180 degrees for the R–HS angle, are separated by a saddle point. The coordinate variables in this potential are RCM , the distance from the rare gas to the center of mass of SH, and θ, the angle between the RCM and rSH vectors, with θ = 0 corresponding to the linear R–HS geometry. The HS/D bond distance, rSH , is fixed at its v = 0 level value in the isolated moiety for this model. The potential following Bowman’s prescription is given by V (RCM , θ) = V0 (RCM )[1 − f (θ)] + Vsp (RCM )f (θ), 0 ≤ θ ≤ θsp

(1)

V (RCM , θ) = Vsp (RCM )[1 − g(θ)] + Vπ (RCM )g(θ), θsp ≤ θ ≤ π.

(2)

where V0 , Vsp , and Vπ , are radial cuts of the potential for θ equal to zero, the saddle point angle, and π respectively. These radial potentials each have a generalized Morse functional form 0

V (RCM ) = D[e−α2 (RCM −Re ) − 2e−α1 (RCM −Re ) ]. 4

(3)

The minimum of this potential is at Re . The value of V (Re ) is equal to −De , where De (the dissociation energy) is the well depth at Re . Hence 1 2α1 ln . α2 α2

0

Re = Re +

(4)

and D=

De . 2(1 − α1 /α2 )

(5)

The functions f (θ) and g(θ) are switching functions that vary between 0 and 1 for the ranges 0 ≤ θ ≤ θsp and θsp ≤ θ ≤ π, respectively. In terms of the variable x, f (x) and g(x) have the same expression, i.e. a b 3a 3a b a f (x) = x2 + (10 + − )x3 + ( − b − 15)x4 + (6 + − )x5 , 0 ≤ x ≤ 1 2 2 2 2 2 2

(6)

where a = fxx (0) and b = fxx (1). In general, x=

eλ1 θ − 1 eλ1 θsp − 1

f or 0 ≤ θ ≤ θsp

(7)

and x=

eλ2 (θ−θsp ) − 1 eλ2 (π−θsp ) − 1

f or θsp ≤ θ ≤ π

(8)

where λ1 and λ2 are free parameters determined in the potential optimization.

B. The Sequential Diagonalization Truncation Method

In this section we discuss the implementation of the IRLM approach by Sorensen et al.13 along with our implementation of the SDT approach by Choi and Light.18 In our previous paper,11 we discussed only the implementation of the IRLM method. Combining these two methods allows us to reduce greatly the amount of computational time. In the original implementation of Choi and Light, the SDT transformation was carried out explicitly by transforming from the DVR representation of the Hamiltonian to the lower order SDT 5

representation. This transformation converts the sparse DVR Hamiltonian into a dense SDT Hamiltonian.18 Since the IRLM approach only requires the action of the SDT Hamiltonian on a vector, the sparsity of the original DVR matrix can be preserved while obtaining the eigenpairs of interest. At this point we need to introduce some additional equations in order to explain the overall efficiency of method. The matrix elements of the full Hamiltonian from Eq.1 of Korambath et al.11 can be written as 0

0

0

αβ K = HαβK

XX j0 j

=

i0 i

dR α0 ,α +

0 0

(Kθ)

0

(Kθ)

i j K R Tj 0 β 0 TiR0 α0 · HijK · Tiα Tjβ

· δβ 0 β · δK 0 K

~2 + 2µ



1 1 + 2 2 Rα re

 (θK) dβ 0 β · δα0 α · δK 0 K

~2 {[J(J + 1) − 2K 2 ] · δα0 α · δβ 0 β · δK 0 K 2 2µRα

+ 0 0 − (1 + δK0 )1/2 Λ+ JK Bβ 0 β · δα α · δK K+1 − 0 0 0 0 0 − (1 + δK 0 0 )1/2 Λ− JK Bβ 0 β · δα α · δK K−1 } + V (Rα , θKβ ) · δα α · δβ β · δK K

(9)

where the V (Rα , θkβ ) matrix is from the potential energy term, the dR matrix is from the radial kinetic energy, the d(θK) matrix is from the angular kinetic energy, the B + and B − matrices are from the coriolis terms. α and β are the indices for the DVR quadrature nodes for radial and angular kinetic energy matrices respectively. The value of K ranges from 0 to J. We carry out the diagonalization of the above Hamiltonian in three stages. Following Choi and Light, we choose the initial diagonalization for the 1D Hamiltonian with respect to the radial coordinate. 1D

0

0

0

αβ K HαβK = dR + V (Rα , θkβ ) · δα0 α · δβ 0 β · δK 0 K + α0 ,α

≡ in which

1D (βK) hα0 ,α

1D (βK) hα0 ,α

· δβ 0 β · δK 0 K

~2 [J(J + 1) − 2K 2 ]δα0 α · δβ 0 β · δK 0 K 2µRα2 (10)

0

are the α , α element of the 1D Hamiltonian matrix in any of the (β, K)

diagonal blocks of the full

3D

H matrix. Each pair of (β, K) defines the values for θβ,K ,

corresponding to the β th DVR point in the angular basis, and K for the projection of 6

total angular momentum J onto the body fixed axis. The diagonal elements V (Rα , θkβ ) and ~2

2µR2α

[J(J + 1) − 2K 2 ] are also included for convenience.

Each

1D (βK)

matrix is then diagonalized using a standard EISPACK or LAPACK19

h

diagonalizer to give 1D (βK)

h

in which

1D

=

1D

C (βK) ·1D E (βK) · [1D C (βK) ]T

C (βK) is a NR × NR 1D-eigenvector matrix, and

1D

(11) E (βK) is a diagonal matrix

containing eigenvalues of the 1D Hamiltonian for each (β, K) blocks. In the SDT approach during this operation, the high energy eigenvectors in

1D

C (βK) were truncated to speed up

the diagonalization either by retaining a constant number of vectors or by retaining the states

1D

C (βK) , which satisfy an energy cutoff condition such that 1D

The truncated

1D

1D E (βK) ≤ Ecutoff

C (βK) matrix, denoted by

1D

C˜ (βK) , has a reduced dimension NR × PβK .

Where PβK corresponds to the number of 1D eigenvectors which satisfy the cutoff limit in the corresponding (βK) block. The next step is to calculate the eigenpairs for the 2D Hamiltonian for each K block. In the truncated 1D eigenvector representation, the 2D Hamiltonian, 0

2D l β(K) hlβ

=

X1D α0 α

2D K

h , is given by

0 0 1 1 (θK) (β K ) (βK) (βK) C˜α0 l0 · {1D hα0 ,α · δβ 0 β + ( 2 + 2 )dβ 0 β · δα0 α } · δK 0 K ·1D C˜αl Rα re

(βK)

= 1D El · δl0 l · δβ 0 β · δk0 K X1D (β 0 K 0 ) 1 1 (θK) (βK) + C˜α0 l0 · {( 2 + 2 )dβ 0 β · δα0 α } · δK 0 K ·1D C˜αl Rα re 0

(12)

αα

0

where, l = 1, 2, ......, PβK , l = 1, 2, ........., Pβ 0 K , and β = 1, 2, ..., Nθ . Up to this point the steps in Choi and Light’s18 SDT transformation and our IRLM/SDT approach are the same. In the next step, we need to obtain a certain number of eigenpairs – the exact number depending on the number of energy levels to be determined. In the original STD implementation, the 2D hK Hamiltonian is formed explicitly and a conventional 7

diagonalizer used to determine the eigenpairs. In our approach, we can avoid assembling the 2D K

dense

h

matrix because the IRLM approach only requires that one compute the action

of this matrix on a vector. This procedure is very efficient because we were able to take advantage of the underlying sparsity of the DVR matrix and because we do not have to store the dense

2D K

h

matrix.

Computation of this matrix–vector operation given by (βK)

˜ + C(

y = [1D El

1 1 (θK) ˜ T + )d 0 C ]x Rα2 re2 β β

(13)

is carried out in four steps as follows: u = C˜ T x

v=(

(14)

1 1 (θK) + 2 ) dβ 0 β u 2 Rα re

(15)

w = C˜ v

y=

1D

(βK)

El

(16)

x+w

(17)

At this point, if J = 0 no additional computations are required. However, for J > 0, ˜˜ in the we follow the original approach of Choi and Light to obtain the Hamiltonian, H truncated 2D-eigenvector basis including Coriolis coupling, as 0 ˜˜ m K = X X2D C˜ (K 0 ) · H ˜ l0β 0 K ·2D C˜ (K) H mK lβK lβm l0 β 0 m0 0

0

β0 β

l0 l

(K) = 2D Em · δm0 m +

X X2D 0

0

X1D

(K ) C˜l0 β 0 m0 · {

0

β β

0

(β K C˜αl0

0

)

α

l l

~2 + [−(1 + δK0 )1/2 Λ+ JK Bβ 0 β · δK 0 K+1 2 2µRα

(K) − 1D ˜ (βK) −(1 + δK 0 0 )1/2 Λ+ Cαl } ·2D C˜lβm JK Bβ 0 β · δK 0 K−1 ] ·

(18)

˜ in the representation of the truncated 1D-eigenvectors where the truncated Hamiltonian, H, is given by 0

0

0

˜l β K = H lβK

X1D 0

0

0

0

0

0

(β K ) α β K 1D ˜ (βK) C˜α0 l0 · HαβK · Cαl

αα

8

and 2D

(K) Em · δm0 m =

X X2D 0

β β

0

0

0

(K ) l β (K) 2D ˜ (K) C˜l0 β 0 m0 ·2D hlβ · Clβm

(19)

0

l l

As in the SDT procedure of Choi and Light18 the dimension of the resulting H DV R is ˜˜ of reduced from (Nmax = NR × Nθ × Kmax ) to the successively truncated Hamiltonian H P ˜˜ is diagonalized using a EISPACK routine. size K 2D NK  Nmax . The truncated H

C. Polynomial Preconditioning

The number of iterations taken by the IRLM routine to obtain convergence is an important factor in the diagonalization speed up. So even if we have an efficient matrix–vector product algorithm the speed up is not guaranteed unless we have a favorable eigenvalue spectra. The IRLM procedure takes a smaller number of iterations when the eigenvalues desired are well separated. In this problem, the eigenvalues that we are interested in are tightly clustered at the lower energy region and the eigenvalues that we do not want are well separated at the higher energy region. This is exactly opposite to what is desired to have faster convergence. Fortunately, we can transform this eigenvalue spectrum using a Chebychev polynomial to a problem that will have the same eigenvectors but a much better eigenvalue distribution. Once the IRLM procedure converges to the desired eigenvectors the original eigenvalues can be obtained easily using the Raleigh Ritz quotient. The details of Chebychev polynomial preconditioning is explained in Korambath et al.11 The combination of DVR, SDT and Chebychev preconditioning are the key factors in getting good performance using the IRLM procedure. In Table 1,

we compare the time taken for DVR/IRLM, SDT/DVR, and

SDT/DVR/IRLM approaches to calculated 10 eigenpairs for the Ar·SH complex. Here, we use Lobatto functions for the radial DVR functions and associated Legendre functions as the angular DVR functions as in the previous paper.11 The times, in seconds, are given on

9

calculations performed on an IBM-3CT workstation for J = 0 to J = 2. For SDT/DVR and SDT/DVR/IRLM methods 100 radial functions are truncated to 70 radial functions. The first row gives the time taken for the IRLM approach, and the third row gives the time taken for the approach where we have combined the advantages of IRLM by Sorensen and SDT by Choi and Light.18 As can be seen the SDT/DVR/IRLM method is more than 7 times faster than the original SDT/DVR method of Choi and Light.18 The difference becomes even more significant as the number of radial and angular functions increases.

D. Fitting Strategy

We have carried out the fit to the experimental energy values and rotational constants simultaneously to arrive at the potential energy surface for the R·SH molecules. We have varied 11 parameters that contribute significantly to improving the fit to the experimental values.

These parameters correspond to the Morse parameters for

0 degrees and saddle point (sp) of the potential energy surface, including specifically Re (0), Re (sp), α1 (0), α1(sp), α2 (0), α2 (sp) in Eq. (3) and De (0) and De (sp) in Eq. (5), and λ1 in Eq. (7) and fxx (0). We also varied the angle for the saddle point. The corresponding parameters for the 180 degree angle were approximated qualitatively to resemble the potential energy surface of Ar·OH where some ab initio information was available for that region of the PES. The starting estimates for the parameters were approximately the same as their optimized values for the Ar·OH system, with appropriate adjustments for the van der Waals radii. The van der Waals stretch (vs ) eigenvalues, arising from the J = 0 calculation and one quanta (v1b ) of bend eigenvalues arising from the J = 1 calculation, were fit to their corresponding experimental values. The rotational constants were calculated from the eigenfunctions as expectation values of < 1/R2 > to compare with the experimental rotational constants. Rotational constants can also be calculated from the difference between eigenvalues of the lowest rotational levels. The results were equivalent within ±0.0001cm−1 for 10

the R·SH complexes. For these calculations the diatom rotational constant Bo = Be − 0.5αe is from Huber and Herzberg.20 Using the optimized surface for the R·SH complexes we have predicted the vibronic energy levels and rotational constants for the R·SD complexes. To accomplish this, we express the mass weighted RCM and θ coordinates for the R·SD complex in terms of those used for R·SH, which are now denoted R0CM and θ0 . This correction is necessary because the center of mass of the diatom shifts towards deuterium by an amount ∆= 0.0408 ˚ A. Hence the RCM and θ coordinates have to be properly shifted when the potential energy is calculated for the R·SD molecule. The new RCM is given by 0 2 0 0 0 (RCM + 2RCM cosθ0 ∆ + ∆2 )1/2 and the new cosθ by (RCM cosθ0 + ∆)/RCM , where ∆ is

the center of mass shift between SH and SD. This procedure was helpful in assigning the vibrational quantum numbers for the deuterium species as there are less experimental data and the experimental absolute vibrational quantum number assignment was not as certain as it was for the hydrogenated species.

III. RESULTS AND DISCUSSION

The results of the nonlinear fit, using the experimental2,3 values to arrive at a global potential, are reported here. For the Ar·SH system we varied the 11 parameters that contributed significantly in improving the fit to 15 experimental values: 7 vibrational energy levels and 8 rotational constants. For the Kr·SH system we again varied the 11 parameters to fit 17 experimental values: 14 vibrational energy levels and 3 rotational constants. Generally speaking, for experimental energy levels, the agreement was usually within ± 0.5 cm−1 , which is within the accuracy of the lower resolution (though not the high resolution ones) and more abundant experimental survey scans, and within ± 0.0001 cm−1 accuracy for the rotational constants. In the previous section we explained the procedure for using the R·SH potential surface to calculate the eigenvalues for the deuterated complex. This approximation works very well for Ar·SH/D, but less well for Kr·SH/D. As in the case of the

11

Ar·OH/D system an inverse isotope effect is observed in both the Ar·SH/D and Kr·SH/D systems. This inverse isotope effect is reproduced by our calculations.

A. Ar·SH/D

The experimental vibrational assignment for the hydrogenated complex was done using the normal isotopic shift relationships of the heavier

34

S and

32

S isotopic species. The

details of this analysis are given in a separate publication. 2 In Table 2, the experimental and calculated values used in the fit for the Ar·SH complex are compared with the vibrationless e state taken as the zero of energy. The eigenvalues all agree within 0.3 cm−1 level of the A with the experimental vibrational energy levels and the calculated rotational constants are generally within 0.0001 cm−1 of the corresponding experimental values. This is a very good agreement considering the nature of the PES used in the optimization process. Table 3, gives all the eigenvalues and rotational constants predicted by the optimized Ar·SH surface up to the (0,0,6) level including the corresponding values for Ar·SD. The average angle for the rotation of H/D atom off-axis, as calculated from the expectation value of < cosθ >, is also given in the last column of the table. In Table 4 eigenvalues intervals (En − En−1 ) and rotational constants predicted for Ar·SD using the optimized Ar·SH surface are compared with the experimental values. Focusing on the difference in the adjacent eigenvalues is more appropriate because the lowest experimental eigenvalue measured for Ar·SD is (0,0,2) and extrapolating back to (0,0,0) level may cause some inaccuracy. In Table 5, all the parameters which determine the PES for Ar·SH are given. As noted in Sec. II D, only the parameters corresponding to 0 degrees and the saddle point were optimized. Fig. 1 shows the potential energy surface for the Ar·SH molecule. This plot has a minimum of -877.2 cm−1 at 0 degrees, a minimum of -202.9 cm−1 corresponding to 92.8 degrees which is called the saddle point, and a minimum of -700 cm−1 (unoptimized) corresponding to 180 degrees. The qualitative nature of the surface resembles other van der Waals systems with minima at 0 and 180 degrees and a saddle point close to 90 degrees. The fit value of 12

the Re at 0 degrees is 6.457 au (3.417 ˚ A), at the saddle point 8.090 au (4.281 ˚ A), and at 180 degrees it is an unoptimized 4.2 au (2.2 ˚ A). The contours, shown with narrow solid lines, are each separated by 50 cm−1 and drawn between -25 and -875 cm−1 . The thick solid line is the zero of energy for in the surface and the dashed lines extend between 100 and 500 cm−1 . In Fig. 2, Fig. 3, and Fig. 5 the probability density corresponding to the (0,0,0), (0,0,5), and (0,1,1) levels respectively, of the hydrogenated complex, are overlapped on the potential energy surface to illustrate the extent each of the eigenfunctions samples the PES. It may be observed that the majority of the PES sampled by the pure van der Waals stretch is in the region ≤ 90 degrees of motion of the inert gas around the SH axis, with the linear configuration R·SH as 0 degrees. Fig. 4 shows the probability density for the (0,0,5) level of the deuterated complex. This can be directly compared with Fig 3, which is the same level for the hydrogenated species. As can be seen from this plot the area of the PES sampled by the deuterated complex does not extend nearly as far away from the linear Ar–D–S geometry as does the hydrogen complex. The wavefunctions for levels with bending motion excited extend much further out into the angular part of the potential, as is shown in Fig. 5. This PES would not be expected to give accurate quantitative information for geometries corresponding to Ar·SH angles near 180 degrees because none of the wavefunctions corresponding to the experimentally observed states have appreciable amplitudes for that region of the potential. The calculated and experimentally approximated D0 values for Ar·SH are given in Table 6. The calculated D0 value is the lowest eigenvalue for the J = 0 calculation and the calculated De value is the minimum of the potential surface. The difference in the experimental and calculated D0 value is due to the fact that all the experimental values are measured far below the dissociation limit. An accurate experimental determination is very difficult in the absence of energy levels close to dissociation limit. In Table 7, the difference in average RCM values, defined as the distance from the inert gas to the center of mass of the SH/D moiety, corresponding to various stretch levels for the

13

Ar·SH/D isotopes are given. Average RCM values were calculated from the rotational constant, using the (< 1/R2 >)−1/2 inversion. The experimentally observed differences, of about 0.1 ˚ A, are nearly 100 times larger than in tight chemically bound molecules. Nonetheless the calculation does an excellent job of reproducing the experimentally observed differences. The calculated values are in good agreement with the two experimentally measured values.

B. Kr·SH/D

The approach taken for Kr·SH is the same as for Ar·SH. Using experimental data for the 86

Kr and 84 Kr isotopes we were able to unambiguously assign the absolute vibrational quan-

tum numbers for the hydrogenated complex, the details of which are presented elsewhere.2 For Kr·SH the lowest experimentally observed level corresponds to the (0, 00 , 5). Rather than using a long extrapolation to the vibrationless level, which could introduce considerable error in the relative frequencies from the true origin, we carried out the fit of the energy surface by minimizing the difference with respect to that of the (0, 00 , 5) level. The vibrational energy values and rotational constants used in the fit and their optimized calculated values are compared in Table 8. The calculated eigenvalues are within ± 0.3 cm−1 (except for two perturbed levels) of the experimentally observed vibrational energy values and the rotational constants are within 0.0001 cm −1 (again except for the perturbed levels). It should be noted that the agreement between the calculated and experimental values given in Table 8 are for frequency differences relative to the (0,00 ,5) level. If we use the experimentally extrapolated origin the error is approximately a constant 6 cm−1 for all the levels. We attribute this to the error incurred by the long extrapolation. Table 9 contains the eigenvalues and the rotational constants calculated from the optimized PES for Kr·SH along with corresponding predicted values for Kr·SD. In Table 10, the calculated rovibrational values (En −E(n−1) ) for Kr·SD as predicted from the optimized Kr·SH surface are given. The lowest experimental value observed for Kr·SD is the level (0,00 ,7). This comparison avoids any error in the extrapolation to the zeroth 14

level. The agreement is not at all bad, but not as good as for Ar·SD prediction. There are several possible reasons for the poorer agreement: (i) the v = 0 adiabatic approximation for the SH and SD diatomic molecules may not hold as well at higher stretch values (towards the dissociation limit) – in this case the first energy levels of the van der Waals complex are about midway between the v = 0 and v = 1 levels of the SD diatomic molecule; (ii) limited flexibility of the potential form in the angular region; (iii) limitation in the number of experimental levels available to carry out the fit; (iv) the assumption that this potential is transferable from SH to SD may not be as good for the higher vibrational levels because this approximation does not account for the differences in the average polarizability and dipole moment of the v = 0 state of the two diatomic molecules. Table 11 lists the optimized parameters for the Kr·SH fit corresponding to the results given in Table 8. Again only those parameters corresponding to 0 degree and the saddle point were allowed to vary. Fig. 6 shows the potential energy surface for the Kr·SH system. The surface has a minima of -1706.2 cm−1 corresponding to 0 degrees, an unoptimized -1400 cm−1 corresponding to 180 degrees. The saddle point for this surface is located at 64.0 degrees where the minimum is -158.9 cm−1 . The Re values are 6.168 au (3.264 ˚ A), 7.589 au (4.016 ˚ A), and 4.2 au (2.2 ˚ A) for 0 degrees, the saddle point, and 180 degrees (unoptimized) respectively. The contours, shown with narrow solid lines, are drawn between -50 and -1700 cm −1 with intervals of 100 cm−1 . The thick solid line corresponds to 0 cm−1 and the dashed lines are drawn between 500 and 2000 cm−1 . In order to illustrate the region of the potential that the experimentally observed levels sample we have plotted the probability density on top of the potential energy surface for some representative levels. These are shown in Fig. 7, Fig. 8, Fig. 9, and Fig. 10. The figures show that most of the wavefunctions sample the surface well up to about 90 degrees from the linear R·SH configuration, and to a limited extent beyond, eg. see Fig. 10. Figure 11 is a plot of the (0,2,6) level of the deuterated complex. This figures depicts the same level shown for the hydrogenated complex in Fig. 10. Here, as with the Ar complex, we see that the area of the potential sampled by the deuterated complex is significantly less 15

than that sampled by the hydrogenated complex. The calculated and experimental D0 values and calculated De value for this molecule are given in Table 6. The calculated D0 value is the lowest eigenvalue for the J = 0 calculation and the calculated De value is the minimum in the potential plot. The table shows good agreement between the calculated values and experimental estimates for D0 . One can note the much improved agreement between calculated and experimental values for Kr·SH vs. Ar·SH where, for the former, levels much closer to the dissociation limit were measured. Table 7, gives the difference in the average RCM values for H and D isotopomers. For the single overlapping stretch level for which we have experimental data, there is good agreement in the bond length difference between the hydrogenated and deuterated species.

IV. CONCLUSIONS

Recent experimental data for the Ar·SH and Kr·SH van der Waals complexes offer the opportunity for the development of meaningful PES. We have used this data to optimize the parameters in an analytical potential model. This approach involves combining the IRLM and SDT methods to take advantage of the sparse nature of the DVR Hamiltonian. This procedure resulted in a significant reduction in the computational time needed to optimize the parameters to fit the data. The area of the potential near 0 degrees, corresponding to a linear R·SH geometry, has been calculated with very good agreement compared to the experimental results. The saddle point for the complexes has also been fit using the experimental data. The area of the potential near 180 degrees has only been qualitatively produced as the experiment does not sample well this region of the PES. Vibrational eigenvalues and rotational constants predicted from the PES agree very well with the observed values for both Ar·SH and Kr·SH. Values for the deuterated isotopomers were predicted from the PES optimized for the lightest hydrogen isotope. The agreement with experiment is very good for the Ar·SD

16

complexes. The results of the Kr·SD complex, which involve higher levels of excitation, are somewhat less satisfactory. As with the Ar·OH/D complexes we have observed an inverse isotope effect experimentally. This effect was reproduced quite well in the calculations. ACKNOWLEDGEMENT: The authors gratefully acknowledge support of this work by the National Science Foundation(NSF) via grant 9320909. Additional support was provided by the NSF in the form of Postdoctoral fellowship grants (ASC-9405161 and ASC-9504071). The authors would also like to thank Dr. Anne McCoy for helpful discussions.

17

REFERENCES 1

(a) J. M. Bowman, B. Gazdy, P. Schafer, and M. C. Heaven, J. Phys. Chem. 94, 2226 (1990); 94, 8858 (1990); (b) U. Schnupf, J. M. Bowman, and M. C. Heaven, Chem. Phys. Lett. 189, 487 (1992).

2

M.-C.Yang, C.C.Carter and T.A.Miller, preceeding paper.

3

C. C. Carter and T. A. Miller, preceeding paper.

4

M. -C. Yang and T. A. Miller, to be published.

5

(a)M. C. Heaven, J. Phys. Chem. 97, 8567 (1993); (b)M. C. Heaven, Ann. Rev. Phys. Chem. 43, 283, (1992), and references therein.

6

M. I. Lester, R. A. Loomis, L. C. Giancarlo, M. T. Berry, C. Chakravarty, and D. C. Clary, J. Chem. Phys. 98, 9320 (1993).

7

T. S. Ho, H. Rabitz, S. E. Choi, and M. I. Lester, J. Chem. Phys. 104, 1187 (1996).

8

M.-L.Dubernet and J.M.Hutson, J. Chem. Phys. 99, 7477 (1993).

9

M. Yang and M. H. Alexander, J. Chem. Phys. 103, 3400 (1995).

10

A. D. Esposti and H. -J. Werner, J. Chem. Phys. 93, 3351 (1990).

11

P. P. Korambath, X. T. Wu, and E. F. Hayes, J. Phys. Chem 100, 6116 (1996).

12

(a) J. V. Lill, G. A. Parker, and J. C. Light, Chem. Phys. Lett. 89, 483 (1982); (b) J. C. Light, I. P. Hamilton, and J. V. Lill, J. Chem. Phys. 82, 1400 (1984); (c) J. V. Lill, G. A. Parker, and J. C. Light, J. Chem. Phys. 85, 900 (1986).

13

(a) D. C. Sorensen, SIAM J. Matrix Anal. Appl. 13, 357 (1992). (b) R. Lehoucq, D. C. Sorensen, and P. A. Vu; ARPACK: Fortran subroutines for solving large scale eigenvalue problems, Release 2.1, available from [email protected] in the scalapack directory. In addition anonymous ftp to ftp.netlib.org can be used to retrieve the files. 18

14

a) Z. Ba˘ci´c, R. M. Whitnell, D. Brown, and J. C. Light, Comput. Phys. Commun. 51, 35 (1988); b) R. M. Whitnell and J. C. Light, J. Chem. Phys. 90, 1774 (1989); c) T. J. Park and J. C. Light, J. Chem. Phys. 90, 2593 (1989).

15

P. Pendergast, Z. Darakjian, E. F. Hayes, and D. C. Sorensen, J. Comput. Phys. 113, 201 (1994).

16

R. A. Friesner, J. A. Bentley, M. Menou, and C. Leforestier, J. Chem. Phys. 99, 324 (1993).

17

M. J. Bramley and T. Carrington, J. Chem. Phys. 101, 8494 (1994).

18

S. E. Choi and J. C. Light, J. Chem. Phys. 97, 7031 (1992).

19

The EISPACK/LAPACK diagonalizer can be obtained from the Fortran library routines available from [email protected] in the eispack/lapack directory. In this study we used EISPACK diagonalizer routines.

20

K. P. Huber and G. Herzberg, “ Constants of Diatomic Molecules ”, 1979, p. 590 Van Nostrand Reinhold Company, New York, New York.

19

Figure Captions Figure 1. Ar·SH A 2 Σ+ interaction potential obtained using the optimized parameters given in Table III. (θ = 0o corresponds to Ar–H–S geometry and θ = 180o corresponds to Ar–S–H geometry ) The thick solid contour line corresponds to 0 cm −1 , the narrow solid lines correspond to negative energy contours between –25 cm−1 and –875 cm−1 in increments of 50 cm−1 , and the dashed lines correspond to positive contours between 100 cm−1 and 500 cm−1 in increments of 100 cm−1 . Figure 2. Ar·SH A 2 Σ+ interaction potential overlapped with the probability density for the (0,0,0) energy level to illustrate the extent of overlap for that wavefunction on the PES. We used 100 radial DVR functions (NR = 100) and 80 angular DVR functions (Nθ = 80) to plot the probability density. The probability density is divided into 3 regions denoted by shades of varying intensity to distinguish the different regions. The darkest region corresponds to densities between 0.3858 and 0.1544, the middle shading has a minimum density of 0.0772, with the lightest region being terminated at a density of 0.0001. Figure 3. Ar·SH A 2 Σ+ interaction potential overlapped with the probability density for the (0,0,5) energy level to illustrate the extent of overlap of the eigenfunction of the maximum stretch level observed experimentally. The probability density is divided into 3 regions denoted by shades of varying intensity to distinguish the different regions. The darkest region corresponds to densities between 0.0980 and 0.0392, the middle shading has a minimum density of 0.0195, with the lightest region being terminated at a density of 0.0001. Figure 4. Ar·SH A 2 Σ+ interaction potential overlapped with the probability density for the Ar·SD eigenfunction of the (0,0,5) energy level to illustrate the region of the potential sampled by the deuterated complex compared to that shown in Fig. 3 for the hydrogenated complex.The probability density is divided into 3 regions denoted by shades of varying intensity to distinguish the different regions. The darkest region corresponds to densities between 0.1640 and 0.0984, the middle shading has a minimum density of 0.0327, with the lightest region being terminated at a density of 0.0001. 20

Figure 5. Ar·SH A 2 Σ+ interaction potential overlapped with the probability density for a bending level, (0,1,1), to illustrate the sampling of the PES by a bend levels on the potential surface. The probability density is divided into 3 regions denoted by shades of varying intensity to distinguish the different regions. The darkest region corresponds to densities between 0.0657 and 0.0263, the middle shading has a minimum density of 0.0131, with the lightest region being terminated at a density of 0.0001. Figure 6. Kr·SH A 2 Σ+ interaction potential obtained using the optimized parameters given in Table IV. (θ = 0o corresponds to Kr–H–S geometry and θ = 180o corresponds to Kr–S–H geometry ) The thick solid contour line corresponds to 0 cm −1 , the narrow solid lines correspond to negative energy contours between –100 cm−1 and –1700 cm−1 in increments of 100 cm−1 , and the dashed lines correspond to positive contours between 500 cm−1 and 2000 cm−1 in increments of 500 cm−1 . Figure 7. Kr·SH A 2 Σ+ interaction potential overlapped with the probability density for the (0,0,0) energy level to illustrate the extend of overlap for that wavefunction on the PES. We used 150 radial DVR functions (NR = 150) and 80 angular DVR functions (Nθ = 80) to plot the probability density. The probability density is divided into 3 regions between 0.4962 and 0.0001 shown by shades of varying intensity to distinguish the adjacent regions. The darkest region corresponds to the region between density 0.4962 and 0.1985, each subsequent region has its density reduced by 0.0992, with the lightest region being terminated at a density of 0.0001. Figure 8. Kr·SH A 2 Σ+ interaction potential overlapped with the probability density for the (0,0,13) energy level to illustrate the extent of overlap of the eigenfunction of the maximum stretch level observed experimentally. The probability density is divided into 3 regions denoted by shades of varying intensity to distinguish the different regions. The darkest region corresponds to densities between 0.0759 and 0.0304, the middle shading has a minimum density of 0.0151, with the lightest region being terminated at a density of 0.0001. Figure 9. Kr·SH A 2 Σ+ interaction potential overlapped with the probability density 21

for a one quantum bend level, (0,1,9), to illustrate the extent of the eigenfunction on the PES. The probability density is divided into 2 regions denoted by the different shades. The dark region corresponds to densities between 0.0519 and 0.0104 with the light region being terminated at a density of 0.0001. Figure 10. Kr·SH A 2 Σ+ interaction potential overlapped with the probability density for the (0,2,6) energy level to illustrate the extent of the overlap of the highest bend level observed experimentally. Note the node in the eigenfunction in the θ= 20-30o region. The probability density is divided into 3 regions denoted by shades of varying intensity to distinguish the different regions. The darkest region corresponds to densities between 0.1065 and 0.0426, the middle shading has a minimum density of 0.0212, with the lightest region being terminated at a density of 0.0001. Figure 11. Kr·SH A 2 Σ+ interaction potential overlapped with the probability density for the Kr·SD eigenfunction of the (0,2,6) energy level to illustrate the difference in the extent of the region sampled by the deuterated eigenfunction compared with that of the hydrogenated complex for the same level (shown in Fig. 10). The probability density is divided into 3 regions denoted by shades of varying intensity to distinguish the different regions. The darkest region corresponds to densities between 0.1826 and 0.0731, the middle shading has a minimum density of 0.0365, with the lightest region being terminated at a density of 0.0001.

22

Table 1. Comparison of run times in seconds for DVR/IRLM, SDT/DVR and SDT/DVR/IRLM methods. Calculation times are for 10 eigenpairs of the rovibrational energy levels of the Ar·SH van der Waals molecule for different J values on an IBM workstation.

Method DVR/IRLM SDT/DVR SDT/DVR/IRLM

J=0 size time 2000 44.6 1400 49.5 1400 6.3

23

J=1 size time 4000 75.9 2800 99.0 2800 13.5

J=2 size time 6000 124.4 4200 148.5 4200 21.1

Table 2. Comparison of experimental vibrational eigenvalues and rotational constantsa with calculated values obtained from the A˜ 2 Σ+ Ar·SH PES.

Level (0,0,0 ) (0,00 ,1) (0,00 ,2) (0,00 ,3) (0,00 ,4) (0,11 ,1) (0,00 ,5) (0,11 ,2)

(En − E0,0,0 ) Expt. Calc. Eexp − Ecal 0.0(5) 0.0 0.0 79.7(5) 79.7 0.0 152.1(5) 151.8 0.3 216.6(5) 216.3 0.3 273.4(5) 273.4 0.0 278.9(5) 278.9 0.0 322.7(5) 322.7 0.0 330.6(5) 330.3 0.3

Bv Expt. 0.077151(100)b 0.074484(10) 0.071527(6) 0.068492(6) 0.065262(7) 0.070857(24) 0.061717(8) 0.066630(20)

All values in cm−1 . b Extrapolated value. a

24

Calc. Bexp − Bcal 0.077128 0.000023 0.074407 0.000077 0.071544 -0.000017 0.068509 -0.000017 0.065246 0.000016 0.070892 -0.000035 0.061628 0.000089 0.066553 0.000077

Table 3. Predicted vibrational energies and rotational constantsa for A˜ 2 Σ+ Ar·SH/D from the calculated PESb .

Level (0,00 ,0) (0,00 ,1) (0,00 ,2) (0,00 ,3) (0,11 ,0) (0,00 ,4) (0,11 ,1) (0,00 ,5) (0,11 ,2) (0,20 ,0) (0,00 ,6) (0,11 ,3)

Ar·SH Energy Bv c 0.0 0.07713 79.7 0.07441 151.8 0.07154 216.3 0.06854 217.6 0.07458 273.3 0.06525 278.9 0.07089 322.7 0.06163 330.3 0.06655 355.4 0.06924 363.8 0.05725 370.5 0.06073

d

15.8 16.5 17.4 18.7 24.9 20.7 27.1 23.9 31.1 39.4 33.8 41.5

Ar·SD Energy Bv c 0.0 0.07820 83.5 0.07560 159.6 0.07290 228.7 0.07008 171.1 0.07664 290.7 0.06712 242.7 0.07358 345.8 0.06396 306.1 0.07028 301.8 0.07451 393.9 0.06049 361.3 0.06664

< θ >d 12.9 13.4 14.1 14.9 19.6 16.1 20.7 17.8 22.3 26.7 20.9 24.7

All values given in cm−1 . b The PES was evaluated from Rmin = 2.0 ˚ A to Rmax = 10.0 ˚ A, with 100 radial DVR points, and 20 angular DVR points. c The rotational constant is the expectation value of 1/(2µR2) evaluated from the computed wavefunctions. d The average angle is calculated from < cosθ > using the computed wave functions. a

25

Table 4. Comparison of experimental and predicted rovibrational valuesa for A˜ 2 Σ+ Ar·SD from the optimized PESb for A˜ 2 Σ+ Ar·SH, for J = 0.

Level (0,00 ,1) (0,00,2) (0,00,3) (0,00,4) (0,00,5) (0,00,6)

(En − En−1 ) Bv Expt. Calc. Eexp − Ecal Expt. 83.5 76.2 70.6(10) 69.0 1.6 0.069749(8) 61.8(1) 62.0 -0.2 0.066682(6) 55.2(10) 55.1 0.1 48.0(10) 48.1 -0.1

Calc. Bexp − Bcal 0.075604 0.072905 0.070085 -0.000336 0.067120 -0.000438 0.063959 0.060492

All values are in cm−1 . b The PES was evaluated from Rmin = 2.0 ˚ A to Rmax = 10.0 ˚ A, with 100 radial DVR points and 20 angular DVR points. a

26

Table 5. Optimized values of the parameters for the Ar·SH/D PES.

Parameter (θ = 0) Re (au) 6.4575 α1 0.7117 α2 2.4947 −1 De (cm ) 877.2 fxx 7.254 gxx λ1 λ2 Saddle point a

(θ = sp) (θ = 180) 8.0911 4.18a 0.2730 1.5a 1.2724 0.7625a 202.9 700.0a a -1.0 6.75a -7.11a -0.6440 1.0000a 92.8

Parameters which are not varied in the fit.

27

Table 6. Comparison of calculated and experimental D0 values, in cm−1 , for Ar·SH/D and Kr·SH/D. Molecule De (calc) D0 (calc) D0 (Expt)a Ar·SH 877.2 572.3 458(1) Ar·SD 877.2 637.3 528(5) Kr·SH 1706.2 1291.6 1273(3) Kr·SD 1706.2 1388.7 1425(22) a

Errors are representative of statistical errors only and do not include limitations in the experimental extrapolation.

28

a Table 7. Experimental and calculated differences for the average RCM values in ˚ A for Ar·SH/D and Kr·SH/D.

Ar·SH/D Kr·SH/D Level ∆RCM (expt.) ∆RCM (calc.) ∆RCM (expt.) ∆RCM (calc.) (0,0,3) 0.0772(3) 0.0714 (0,0,4) 0.0847(3) 0.0835 (0,0,11) 0.0988(9) 0.0946 a

Average RCM value is calculated as (< 1/R2 >)−1/2 from the computed wavefunction.

29

Table 8. Comparison of experimental vibrational eigenvalues and rotational constantsa with calculated values obtained from the A˜ 2 Σ+ Kr·SH PES.

Level (0,00 ,6) (0,00,7) (0,20,2) (0,00,8) (0,00,9) (0,11,6) (0,00,10) (0,11,7) (0,20,4) (0,00,11) (0,11,8) (0,20,5) (0,00,12) (0,11,9) (0,20,6) (0,00,13)

(En − E0,0,5 ) Expt. Calc. Eexp − Ecal 94.8(10) 94.9 -0.1 182.6(10) 182.6 0.0 b 260.2(10) 267.0 -6.8 264.6(10)b 263.3 1.3 336.8(10) 336.8 0.0 347.0(10) 347.3 -0.3 403.5(10) 403.4 0.1 419.6(10) 419.3 0.3 431.1(10) 431.0 0.1 463.1(10) 462.9 0.2 483.5(10) 483.3 0.2 498.6(10) 498.7 -0.1 515.7(10) 515.5 0.2 538.7(10) 538.9 -0.2 554.4(10) 554.3 0.1 561.3(10) 561.1 0.2

All values in cm−1 . Energy level is perturbed.

a b

30

Bv Expt.

Calc.

Bexp − Bcal

0.057950(25) 0.058423 0.054459(9) 0.053204 0.051330(21) 0.051434

-0.000473 0.001255 0.000104

0.049484(17) 0.049562

-0.000078

0.047471(20) 0.047562

-0.000091

Table 9. Prediction for vibrational energies and rotational constantsa for A˜ 2 Σ+ Kr·SH/D from the calculated PESb . Level (0,00 ,0) (0,00,1) (0,00,2) (0,11,0) (0,00,3) (0,11,1) (0,00,4) (0,11,2) (0,00,5) (0,20,0) (0,00,6) (0,20,1) (0,11,4) (0,00,7) (0,00,8) (0,20,2) (0,11,5) (0,00,9) (0,11,6) (0,20,3) (0,00,10) (0,11,7) (0,20,4) (0,00,11) (0,11,8) (0,20,5) (0,00,12) (0,11,9) (0,20,6) (0,00,13)

Kr·SH Energy Bv c 0.0 0.06492 132.1 0.06363 256.5 0.06230 336.1 0.06345 373.5 0.06093 455.1 0.06205 483.1 0.05951 566.2 0.06059 585.3 0.05803 653.2 0.06170 680.2 0.05649 757.0 0.06012 764.9 0.05746 767.9 0.05489 848.5 0.05320 852.2 0.05842 852.7 0.05577 922.1 0.05143 932.5 0.05396 938.8 0.05659 988.6 0.04956 1004.6 0.05201 1016.3 0.05455 1048.2 0.04756 1068.5 0.04987 1084.0 0.05216 1100.8 0.04539 1124.1 0.04743 1139.6 0.04869 1146.4 0.04299

d

12.7 13.1 13.5 18.3 13.9 19.0 14.4 19.6 14.9 23.2 15.5 24.2 21.2 16.1 16.8 25.4 22.2 17.7 23.4 26.9 18.8 24.8 28.9 30.5 26.9 32.2 21.9 30.1 43.3 24.7

Kr·SD Energy Bv c 0.0 0.06551 134.1 0.06425 260.8 0.06295 243.5 0.06448 380.2 0.06161 368.5 0.06315 492.5 0.06023 486.0 0.06177 597.6 0.05880 482.1 0.06332 695.7 0.05731 597.2 0.06190 698.9 0.05886 786.9 0.05577 871.3 0.05417 704.6 0.06043 794.4 0.05732 948.9 0.05250 882.6 0.05570 804.4 0.05887 1019.8 0.05075 963.7 0.05400 896.5 0.05724 1084.1 0.04891 1037.6 0.05220 981.0 0.05551 1141.9 0.04697 1104.3 0.05028 1057.7 0.05366 1193.2 0.04488

< θ >d 10.8 11.1 11.5 15.5 11.8 15.9 12.2 16.5 12.7 19.2 13.1 19.8 17.7 13.6 14.2 20.6 18.3 14.8 19.1 21.4 15.6 20.0 22.3 16.4 21.0 23.4 17.5 22.2 24.7 18.9

All values are in cm−1 . b The PES was evaluated from Rmin = 2.0 ˚ A to Rmax = 10.0 ˚ A, with 160 radial DVR points, and 30 angular DVR points. c The rotational constant is the expectation value of 1/(2µR2) evaluated from the computed wavefunctions. d The average angle is calculated from < cosθ > using the computed wave functions. a

31

Table 10. Comparison of experimental and predicted rovibrational valuesa for A˜ 2 Σ+ Kr·SD from the optimized PESb for A˜ 2 Σ+ Kr·SH for J = 0.

Level (0,00 ,8) (0,00,9) (0,00,10) (0,00,11) (0,00,12) (0,00,13)

(En − En−1 ) Bv Expt. Calc. Eexp − Ecal Expt. Calc. Bexp − Bcal 87.5(10) 84.4 3.1 0.054172 81.0(10) 77.6 3.4 0.052500 73.2(10) 70.9 2.3 0.050751 66.7(10) 64.3 2.4 0.048663(44) 0.048914 -0.000251 59.5(1) 57.8 1.7 0.047285(52) 0.046969 0.000316 54.5(10) 51.3 3.1 0.044880

All values are in cm−1 . b The PES was evaluated from Rmin = 2.0 ˚ A to Rmax = 10.0 ˚ A, with 160 radial DVR points and 30 angular DVR points. a

32

Table 11. Optimized values of the parameters for the Kr·SH/D PES.

Parameter (θ = 0) Re (au) 6.1686 α1 0.9581 α2 2.7668 −1 De (cm ) 1706.2 fxx 7.4308 gxx λ1 λ2 Saddle point a

(θ = sp) (θ = 180) 7.5892 4.14a 0.5212 1.85a 2.6900 0.8625a 158.9 1400.0a -1.0a 6.75a -7.11a 0.8269 1.0a 64.0

Parameters which are not varied in the fit.

33