Methods for plotting trajectories for triatomic molecules and visualizing

Jul 1, 1983 - D. W. Noid, D. M. Wardlaw, M. L. Koszykowski, R. A. Marcus. J. Phys. Chem. , 1983, 87 (15), pp 2733–2738. DOI: 10.1021/j100238a011...
1 downloads 0 Views 1MB Size
J. Phys. Chem. 1903, 87,2733-2738

coupled. If, for example, coriolis coupling between vibrational and rotational states is extremely weak, which might be true for small J , then the Af/Af’ symmetry of vibrational and rotational states would be conserved separately, not just the composite ro-vibrational symmetry. In the tunneling region of reaction 1.1,for example, there would then be a factor of 20 between rate constants for states whose vibrational symmetries are A’ and A”. Symmetry alone, of course, cannot answer the question of whether or not there exist weak coupling, i.e., weak selection rules, between states within a given irreducible

-

2733

representation; for this one needs a Hamiltonian and dynamics.

Acknowledgment. This work has been supported by the Director, Office of Energy Research, Office of Basic Energy Sciences, Chemical Sciences Division of the US.Department of Energy under Contract DE-AC03-76SF00098. All calculations were carried out on a Harris H800 minicomputer supported by the National Science Foundation, Grant CHE-79-20181. Registry No. H2C0, 50-00-0.

Methods for Plotting Trajectories for Triatomic Molecules and Visualizing the Caustics D. W. Noid,+ D. M. Wardlaw,* M. L. Koszykowskl,i and R. A. Marcus” Oak RMge National Laboratory, Oak RMge, Tennessee 37830, Arthur Amos Noyes Laboratory of Chemical Physics, California Institute of Technology. Pasadena, Californ& 9 1125, and Sandia National Laboratory, Livermore, California 94550 (Received: January 3, 1983)

A simple and efficient computational method is introduced for visualizing the structure of three-dimensional classical trajectories. Examples are presented for Hamiltonians with incommensurate frequencies, as well as with 1:l and 2:l resonances. Results are also presented for several states of OCS, one of which shows an accidental 41 resonance. Implicationsof this method for path integral quantization and perturbation methods me discussed.

Preface Henry Eyring was one of the pioneers in the recent theory of chemical reaction rates. R.A.M. was fortunate to know him for some three decades, and to see the joyousness with which he approached his science and his interactions with people. In recent years there has been an increasing amount of interest in the field of intramolecular dynamics and the anharmonic motion of vibrations, subjects related to one of Henry’s interests-the statistical theory of mass spectral decomposition patterns. It is a pleasure to dedicate this article to his memory. Introduction Nonlinear dynamics has been the subject of intense study in recent years, and semiclassical methods have been developed to treat bound and quasibound states of multidimensional nonseparable systems.’J Most studies have focused on coupled two-oscillator systems. The two principal methods used for describing the trajectories in two-oscillator systems have been (1)plots of the trajectories themselves in various coordinates and (2) plots of Poincar6 surfaces of section in various coordinates. It appeared at first that both of these techniques would be difficult to extend to three-oscillator systems. The surface of section + Oak

Ridge National Laboratory. Authur Amos Noyes Laboratory of Chemical Physics. NSERC of Canada Postdoctoral Research Fellow. Sandia National Laboratory. *Author to whom correspondence should be addressed at the Author Amos Noyes Laboratory of Chemical Physics. Contribution No. 6769. f

Q022-365418312Q87-2733$Ql.SO10

technique has been extended to three-coordinate systems for maps3 and for trajectorie~,~ and can similarly be extended to four coordinates. In the present paper we extend the trajectory plots to three dimensions in a way which permits visualization of shapes for use in semiclassical quantization and for identification of resonances. Examples of three-oscillatorsystems include nonrotating linear and nonrotating nonlinear triatomic molecules, as well as a time-independent form of a two-laser driven o~cillator.~ Although the 3-D Poincar6 surface of section technique4 is useful in studying the dynamics, there is a prerequisite for performing the semiclassical quantization: suitable paths for evaluation of the phase integrals depend on the shape of the trajectories, namely, on the caustics, which are, in turn, affected by internal resonances. It is highly desirable, therefore, to have a meaningful picture of the several-dimensional trajectory itself, in particular one which depicts the caustics. In this exploratory paper a method is devised for this purpose. It is described in section 11, and several applications are given in section 111. The method is expected to be useful for quantization methods which use the trajectories themselves. It should also be useful for perturbation methods in which one tacitly assumes the nature of the caustics, which them(1) D. W. Noid, M. L. Koszykowski, and R. A. Marcus, Annu. Rev. Phys. Chem., 32, 267 (1981). (2) M. Tabor, Adu. Chem. Phys., 46, 73 (1981); P. Brumer, ibid., 47, 201 (1981); S. A. Rice, ibid., 47, 117 (1981). (3) C. Froeschle, Astrophys. Space Sci., 14, 110 (1971). (4) D. W. Noid, M. L. Koszykowski,and R. A. Marcus, J. Chem. Phvs.. _ . 73, 391 (1980).

(5) J. R. Stine and D. W. Noid, J. Phys. Chem., 86, 3733 (1982).

0 1983 American Chemical Society

2734

The Journal of Physical Chemisfry, Vol. 87, No. 15, 1983

Noid et al. 100

selves depend on the presence and spatial structure of resonant trajectories.

Methods By solving Hamilton’s equations Qi = aH/aPi P i = -aH/aQi (1) with subroutine DEROOT? one can calculate, as described below, the values of (Ql, Q2)for which Q3 is some constant value. In this way, one can find data for a plane (Q3 = constant) of a 3-D trajectory in (Q1, Q2, Q3) space. Because the display of these (Ql, Q2)as points would require a large number of such points to visualize the structure of a trajectory, the values of (Ql, Q2) are binned in a two-dimensional array. When a quasiperiodic trajectory is propagated to increasingly longer times, the number of filled bins in this (Ql, Q2)space begins to stabilize, i.e., the trajectory enters no new bins. The technique is readily implemented: DEROOT automatically determines values for which some function G (defined by the user) is zero. For the Q3 = 0 plots, we set G = Q3. For our purposes, we consider 100 bins of Q1 and 100 bins of Q2. In the plots we use one symbol (+) to denote a filled bin in (Q1, Q2) space, and either use another symbol (0) [Figures 1-31, or a contour where the potential energy equals the total energy [Figures 5-81, to indicate the classically forbidden region. In this way, a cross section of the 3-D trajectory can be plotted. This information on which bin is filled for any Q3 is stored in an integer array with (102)2numbers. (Even a minicomputer can be used. If storage is a problem the information can be packed into words.) It is also useful to make a second type of plot, using the approach of Ashton and M~ckerman,~ by constructing (Ql, Q2) cross sections a t different values of Q3. In this case, we set

and solve for the zeroes of this G using the DEROOT program. As before, the information is stored in an integer array, but now with (102)3values when 100 Q3’s are considered. The collection of bins can be plotted in a 3-D perspective by using the plotting package DISPLA~ (we use below a symbol 0 to denote a filled bin). To visualize the trajectory we have found it useful to change the viewpoint by rotating past the object in the Q1Q2plane at a fixed distance and height above it. These several views from different angles will be termed a “fly-by”.

Model Hamiltonians and Results Several model Hamiltonians are used below for illustration, the first being

Trajectories for the case in which wl, w2, and w3 are incommensurate are considered first, with w1 = 0.7, w2 = 1.3, w3 = 1.0, X = -0.1, and 7 = 0.1. The slice Q3 = 0 is plotted in Figure 1for the zeroth order state ( l , O , l ) in ref 4. The trajectory is seen to be “boxlike”, as expected. For the case w1 = 1.4, w2 = 0.7, and w3 = 1.0 the trajectory exhibits a 2:l Fermi resonance in the Q3 = 0 plane. (6) L. F. Shampine and M. K. Gordon, “Computer Solution of Ordinary Differential Equations: The Initial Value Problem”, W. A. Freeman, San Francisco, 1975. (7) C. J. Ashton and J. T. Muckerman, J. Phys. Chem., this issue. (8) The DISPLA plotting package is prepared by Integrated Software Systems Corporation, San Diego, CA 92121.

Y

50-

,h

n ”

50

0

100

X

Figure 1. View of the Q,, Q 2 plane at Q3 = 0 for an incommensurate Hamiltonian. The labels X, Y, and Z correspond to Q,, Q,, and Q3, respectively. I n Figures 1-3 the initial conditions were determined from the unperturbed Hamiltonian as follows: total energy E = E , E , E, with E, = ( n , ’/,) w,, Q, = 0, and Pi = (2Ei)”2,i = 1 , 2 , 3. Here E = 3.2 and ( n , , n 2 , n 3 )= ( l , O , l ) .

+ +

+

50

Y

3

-50 - 5L

0

53

X

Figure 2. Same as Figure 1 , but with a Fermi resonance interaction; E = 2.5, initial actions were not chosen to correspond to integers.

This slice is shown in Figure 2 for arbitrary zeroth order actions and is of a shape similar to those found for 2-D 2:l resonant systemsg (formed from librating figure eights). Thus, the technique shows that the quantization method used in those systems could be applied to the present system. In Figure 3 is shown a Q3 = 0 slice for the case w1 = w2 = w3 = 1. The trajectory, which has arbitrary zeroth order actions, is of the librating type and is similar to those found in the 1:l 2-D system.1° This trajectory also resembles the local mode type found by Lawton and Childll in a (9) D. W. Noid, M. L. Koszykowski, and R. A. Marcus, J. Chem. Phys., 71, 2864 (1979). (10) D. W. Noid and R. A. Marcus, J. Chem. Phys., 67, 559 (1977). (11) R. Lawton and M. S. Child, Mol. Phys., 37, 1799 (1979); R. M. Hedges and W. P. Reinhardt, Chem. Phys. Lett., 91, 241 (1982).

The Journal of Physical Chemisfry, Vol. 87, No. 15, 1983

Methods for Visualizing 3-D Trajectories

0.30

1

I

a

m

0-

0

Y

2735

Q

-0.30 -0.30

0.30

Q2 100

50 X

“0

0.30

I

b

Figure 3. Same as Figure 1, but for a librating trajectory of a 3-D Henon-Heiles Hamiltonian; E = 5.0, initial actions were not chosen to correspond to integers.

f\

I

c

0

0

I -0.30 -c 0

0

30

Q3 0.? C

N

0 0

(b)

-0.

3

0

0.30

Q1

Figure 5. Three views of the (2,2’,3)state of OCS which show a complicated resonance structure. I n Figures 5-8 the initial conditions were determined as described in the caption for Figure 1, except that E, = ( n , l ) ~ , . Here E = 11083.

+

Figure 4. Fly-by views of the trajectory in Figure 3. The successive views correspond to counter-clockwise rotation of the cube about its Z axis.

trajectory can clearly be seen. In Figures 5-7 are plotted trajectories for a model of nonrotating OCS, using the Hamiltonian in ref 12. Each figure provides views of a trajectory in the three planes

two-coordinate system. A 3-D “fly-by” of this trajectory is given in Figure 4, and the spatial orientation of the 3-D

(12) A. Foord, J. G.Smith, and D. H.Whiffen, Mol. Phys., 29, 1685 (1975).

(d)

(e 1

Noid et al.

The Journal of Physical Chemistty, Vol. 87, No. 15, 1983

2738

I

, I

I

-0.33

I 1

02

QZ 0.3: b

I

I

I I

i

0'

Qt

0

I

-0.33

-(

-0.37 - 0 37

I

I

C

i ?'

Q3 0.37 I

3

0

33

Q3

0.33

I

i -0.3;

83

I

i

0

0.33

at -0.37

0

7

Ql

Flgure 6. Three views of a quasibound (5,0°,5) state of O W E = 16848.

defined by Q1 = o, Q2 = o, and Q3 = 0 , and labeled a, b, and c, respectively. The solid line is the equipotential contour determined by the total energy. Figure describes the (2,p,3) zeroth order of OCS. The trajectory displays a higher order, accidental 4:l resonance in the Q2and Q3 degrees of freedom. The 4:l character of this resonance is demonstrated by the double

Flgure 7. Three views of a "chaotic" trajectory of OCS (state 2,4',3) at -180 Q 2 periods; E = 12131.

figure eight pattern of the Q2-Q3 cross section (Figure 5a) and the fact that w3(=1967.0)/w2(=491.7)= 4.00; the frequencies were obtained from a spectral analysis13(Fourier transform) of the coordinates Q1(t), Q2(t),and Q3(t). Although the zeroth order frequencies, w! = 2092.46 and w i = 523.62, are also close to the ratio 4:1, we note that, for (13) D. W. Noid, M. L. Koszykowski, and R. A. Marcus, J. Chem. Phys., 67, 404 (1977).

The Journal of Physical Chemistry, Vol. 87, No. 15, 1983 2737

Methods for Visualizing 3-D Trajectories

most zeroth order OCS states, the addition of the perturbative terms to Ho12removes this zeroth order near degeneracy resulting in incommensurate fundamental frequencies and therefore distorted box-type trajectories. (These results are not presented here.) For the (2,2O,3) state considered here, the initial conditions and the perturbation have combined to restore the 4:l resonance accidentally. The appearance of four “holes” in the Q1-Q3 cross section for this trajectory (Figure 5c), even though the resonance is w3 N 4w2 rather than w3 N 4wl, is explained as follows: Since the Q2-Q3 resonance condition holds for all values of Q1, we expect the 4:l pattern, albeit distorted, to persist over a range of Q1 values. That is, adherence to the condition w3/w2 = 4.0 coupled with the incommensurability of wl(=817.0) results in excluded volumes of configuration space rather than the excluded areas associated with resonances in two degrees of freedom systems. Here the excluded volumes can be pictured by orienting Figure 5a at 90° with respect to Figure 5c and lining up the Q2 axes. Figure 5b simply displays the (solid) Q1-Q3 cross section in between the two innermost excluded volumes. Figure 6 describes a trajectory for the zeroth order (5,0°,5) state of OCS obtained by using the same Hamiltonian. The trajectory, whose energy exceeds that needed for dissociation, is a quasibound one and is of the distorted box type. This trajectory is the 3-D analogue of a resonance state as found in ref 14. In Figure 7 are plotted three views of a bound “chaotic” trajectory of OCS a t 180 Q2 vibrational periods. [The maximum forward time from which this chaotic trajectory could be successfully back integrated15was only 180 Q2 vibrational periods (- 11.7 ps).] The trajectory was detected to be unstable both by its sensitivity to integration error parameters and the associated back integration difficulties, and by its “grassy” power spectrum.13 We also used the spectral method13 to show that the trajectories of Figures 5 and 6 were quasiperiodic. The trajectory plots in Figure 7 are seen to be somewhat similar to those in Figures 5 and 6, but those in Figure 7 appear to be ragged around the edges. We believe that if executed to sufficiently long times the Figure 7 plots might eventually fill up all the classically allowed space. This point was tested by comparing the time evolution of cross sections for a quasiperiodic (2,0°,3) and a chaotic (2,4O,3)trajectory to a maximum time of 1300 Q2 periods. Despite the onset of numerical integration error, even for the quasiperiodic trajectory, one feature became prominent as time increased. The quasiperiodic cross sections “converged” and retained well-defined boundaries (i.e., the trajectory remained on a torus); a comparison of Figure 8a for 1300 Q2 periods and Figure 8b for 160 Q2 periods illustrates this feature. In the chaotic case the boundaries remained ragged and the trajectory slowly filled the available configuration space (cf. Figure 8c for 1300 Q2 periods and Figure 7c for -180 Q2 periods). The filling of Figure 7a occurred most rapidly, followed by Figure 7c and then Figure 7b.

-

N

0

0-

-0 30 -039

N

0

0 0,

K)

0-

-

033. I

C

-

-

-

-

Conclusions In this paper, a technique is presented for visualizing both the shape and internal structure of classical trajectories. When this technique is coupled with the surface (14) D.W.Noid and M. L. Koszykowski, Chem. Phys. Lett., 73,114 (1980). (15) Back integrated in the chaotic case to at least two digits in the coordinates and one digit in the momenta. (The accuracy in the quasiperiodic cases was considerably greater, as was the forward integration time of -625 Q2vibrational periods in Figures 5 and 6.)

-0331 -033

1

0 Q1

-

-

Figure 8. Comparlson of the quasiperiodic (2,0°,3) state of OCS (a) at 1300 Q 2 periods, (b) at 160 Q 2 periods, and (c) the chaotic (2,4’,3) state at 1300 Q 2 periods; for the (2,0°,3) state E = 10036.

-

of section quantization or with perturbation methods which take the observed shape into account, it is possible to quantize essentially all nonrotating triatomic molecules for which quasiperiodic trajectories are found. The advantage of this technique is that the programming effort is trivial. There will be some cases where the trajectories are “chaotic” even if the quantum mechanical states are “regularly spaced”.16 For such trajectories the present (16) D.W.Noid, M. L. Koszykowski, M. Tabor, and R. A. Marcus, J. Chem. Phys., 72,6169 (1980).

2730

J. Phys.

Chem. 1903, 87,2738-2744

method would, of course, not be useful in the quantization. Another method has been recently developed by Ashton and Muckerman7 which examines the outer surface of the 3-D classical trajectory. The two methods are complementary in that Ashton and Muckerman could obtain our results by “slicing” their representation of a trajectory (in the absence of internal resonances) and we could obtain their results by “stacking” our slices. (The stacking was used to obtain Figure 4.) The resolution of our figures can be altered simply by changing the bin size. In addition, although we have arbitrarily chosen planes perpendicular to the Cartesian axes,

any slice can be obtained by specifying the associated functional form of G in DEROOT. Acknowledgment. We are pleased to acknowledge the support of the present research by the U.S.Department of Energy under Contract W-7405-eng-26with the Union Carbide Corporation (at Oak Ridge) and by a grant from the National Science Foundation (at the California Institute of Technology). D.W.N. is also indebted to Mr. M. A. Bell at Oak Ridge National Laboratory for helpful discussions on the use of the DISPLA plotting package. Registry No. OCS, 463-58-1.

On the Structure of Classical Trajectories in Multldimensional Bound Molecular Systems C. J. Ashton’ and J. T. Muckerman Chemistry Department, Brookhaven National Laboratory, Upton, New York 11973 (Received: February 75, 1983)

A computational method is described which enables visualization of the coordinate space envelopes of classical trajectories in multidimensionalbound molecular systems. The method is exemplified by application to a realistic three-dimensional model of the vibrating water molecule, and its utility in the application of semiclassical quantization techniques is emphasized.

I. Introduction There is considerable current interest in the classical dynamics of nonseparable bound molecular systems,14 and in quasibound and “periodic orbit” trajectories in unbound system^."^ Work in these areas has yielded significant new insights into the physics of nuclear motion on realistic potential surfaces, and semiclassical approximation^^-^*^ enable the extraction from such studies of quantitative information on the quantum mechanics of the systems. Most previous work has been concerned with model systems in two coordinate dimensions: for example, collinear reactive scattering or triatomic molecules with only the stretching vibrations considered. In studying such systems, it has been found essential to visualize graphically the classical motion, both in order to understand its physical significance and to assess the applicability and mode of use of the various semiclassical approximations. In particular, plots of the coordinate space paths of trajectories and of Poincar6 surfaces of section have been found extremely useful. For example, coordinate space plots reveal the shapes of the caustics, which must be known for the application of semiclassical quantization schemes. When systems of three or more dimensions are consid(1)M. Tabor, Adv. Chem. Phys., 46,73 (1980). (2)S.A. Rice, Adv. Chem. Phys., 47,117(1981);P. Brumer, ibid., 47, 201 (1981). (3)I. C. Percival, Adu. Chem. Phys., 36, 1 (1977). (4)D. W. Noid, M. L. Koszykowski, and R. A. Marcus, Annu. Rev. Phys. Chem., 32,267 (1981). ( 5 ) E. Pollak and M. S. Child, Chem. Phys.,60, 23 (1981);E.Pollak and R. E. Wyatt, J. Chem. Phys., 77, 2689 (1982). (6)R. M. Hedges and W. P. Reinhardt, Chem. Phys. Lett., 91,241 (1982). (7)D.W. Noid and M. L. Koszykowski, Chem. Phys. Lett., 73, 114 (1980). (8)R. J. Wolf and W. L. Hase, J. Chem. Phys., 73,3779 (1980). (9)E.Pollak, J. Chem. Phys., 74,5586(1981);Chem. Phys. Lett., 80, 45 (1981).

ered,sJOJ1understanding the classical motion becomes much more difficult. While techniques have been developed for generating surfaces of section in several dimensions,” the problem of visualizing the trajectories themselves has not previously been addressed. The purpose of this paper is to present a practical computational method for visualizing shapes of many-dimensional trajectories. The method is described in section 11, and section I11 presents examples of its use to generate perspective and cross-sectional views of the coordinate space envelopes of three-dimensional trajectories. In section IV we discuss other possible applications and compare the method with another current approach to the problem.

11. Method Method We describe our method as applied to the coordinate space path of a trajectory in a three-dimensional potential well. The first step is to choose a suitable origin in the coordinate space. The choice is not usually critical, any point expected to be in or near the envelope of the trajectory will usually suffice. The angular space around this point is then divided into a large number (n)of cells, each enclosing approximately 4rln steradians. A convenient way to define such a mesh utilizes the conventional polar angles B and 4 relative to the chosen origin and an arbitrarily oriented set of Cartesian axes. The range of 0 (0 to ). is divided into m equal intervals of width r / m , so that interval k ranges from B = ( k - 1 ) r l m to 8 = k r l m ( k = 1, 2, ..., m). In each B interval, the range of 4 (0 to 2.) is divided into Pk intervals of width 2 r / p k , where (10)B.A.Ruf, J. T. Muckerman, C. J. Ashton, and D. W. Noid, to be submitted to J. Chem. Phys. (11)D. W. Noid, M. L. Koszykowski, and R. A. Marcus, J. Chem: Phys., 73,390 (1980).

0022-3654/83/2087-2738$01.50/00 1983 American Chemical Society