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 t~ 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
The Journal of Physical Chemistty, Vol. 87, No. 15, 1983 2739
Classical Trajectories in Multidimensional Molecular Systems (0)
from
\
(c>
from
e=
,
$= 63O
(b)
from
e= 56' ,. 4. = 1 5 3 O
8=124' , 4=243'
(d)
from
8=124'
56'
1 M
, 4333'
/
Figure 1. Perspective views of a quasiperiodic normal-mode trajectory at a total energy of 4804 cm-I. The widely spaced cross wires outline the potential contour at this energy, and the axes labeled 1, 2, and 3 are respectively 6r, (A), 6r, (A), and 6a (units of T ) . The viewing direction for each picture is specific by polar angles (8,4)as defined in the text. See Figure 5 for scaling information.
Pk = max (3, NINT( 2 m sin[ '(
(1)
(and NINT denotes "nearest integer"), so that the cell labeled k j has a range in 8 as described above and a range in 4 from 4 = 27rG - l ) / p k to 4 = 2 r j / p k 6 = 1 , 2 , ...,Pk). The total number of cells is thus
and each cell encloses approximately r 2 /m2 steradians. The time evolution of the selected trajectory is now followed by the usual technique of numerical integration of Hamilton's equations. For each angular cell kj, a record is kept of the largest and smallest distances from the chosen origin, Rkjmax and Rkjmin,respectively, at which the trajectory passes through that cell. After each integration time step, the coordinate space position of the trajectory is examined to determine which angular cell it lies in, and
the maximum and minimum distance records for that cell are updated. If the integration time step is large enough that the trajectory passes through several cells in a given step, it may be desirable to estimate intermediate points on the trajectory path by interpolation. The trajectory is allowed to evolve for some predetermined time, or until the arrays of maximum and minimum distances cease to change significantly. The final result is a mesh of maximum and minimum distance points in coordinate space
k = 1 , 2 ,..., m; j = 1 , 2 ,..., p k
2740
The Journal of Physical Chemistry, Vol. 87, No. 15, 1983
Ashton and Muckerman
(b)
(c)
from
0 = 1 24’ ,. 4. = 2 4 3 ’
t/
M
/
(d)
\
from
e=
,
56’
$=153’
AI
f r o m 0=124’
tI
, 4=333O
M
Figure 2. Perspective views of a quasiperiodic local-mode trajectory at a total energy of 22275 cm-I. See Figure 6 for scaling information.
(where the angular coordinates correspond to the geometrical center of each cell, and the coordinates are relative to the chosen origin and Cartesian axes). These points are then joined together to define a continuous three-dimensional surface (or surfaces) which approximates the envelope of the trajectory. The algorithm used for connecting the points is most easily described in the case where all mesh points are accessed, so that the surfaces of maximum and minimum distance are physically disconnected from each other. Each surface is then constructed separately according to the following conceptual scheme. (a) Each mesh point is joined to its two neighboring points with the same 0 value. (b) Each mesh point is joined to the points on the constant-0 circles above and below it (if these circles exist) which have 4 values closest to that of the point considered. (c) A surface has now been defiied which consists mainly of three-dimensional quadrilaterals, with some triangles. Each quadrilateral is bisected across whichever diagonal subtends the smaller angle at the origin. Thus the surface is finally defined by a set of three-dimensional triangular patches whose vertex coordinates are stored for later use.
When some mesh points are unaccessed, the above algorithm is modified. When the algorithm calls for joining together two unaccessed mesh points, the instruction is ignored. When the algorithm calls for joining an accessed mesh point to an unaccessed one, this is not done but instead the points of maximum and minimum distance for the accessed mesh point are themselves joined together. After bisection of the resulting quadrilaterals, the result is a continuous surface (or surfaces) which is again defined by a set of triangular patches. The resulting surfaces may be examined in several ways. Standard graphical techniques12may be used to produce perspective views. Alternatively, one may generate and plot cross sections in any desired plane, by interpolating line segments in the required plane across all triangles which cross it. 111. Results The utility of our technique is well illustrated by reference to a recent study of classical vibrations of the ro(12) w. M. Newman and R. F. Sproull, “Principles of Interactive Computer Graphics”, McGraw-Hill, New York, 1973.
The Journal of Physical Chemistry, Vol. 87, No. 75, 1983 2741
Classical Trajectories in Multidimensional Molecular Systems
(01
from
\
I (c)
e=
56'
,
+= 63'
(b)
from
AI
f r o m 8=124'
t/
M
e=
56O
,
+=153O
M
,
4=243'
/
(d)
from
8=124'
,. 4=333' .
Figure 3. Perspective views of a quasiperiodic normal-mode trajectory at a total energy of 19 513 cm-'. See Figure 7 for scaling information.
natural coordinate space of the potential. This space is tationless water molecule.1° This system has three coordefined by a pair of O-H bond extensions 6rl and 6r2,and dinate degrees of freedom (one bending and two stretching the H-O-H bending displacement angle 6a. The origin vibrations). The study employed a modified version of the realistic Sorbie-Murre11 potential for the m o l e ~ u l e . ~ ~ ~ lof~ the - ~ ~coordinate system is a t the potential minimum of the molecule, and the sign of 6a is defined such that the Hamilton's equations were formulated and integrated in linear H-O-H geometry corresponds to positive 6a. The the conventional normal coordinates Q1, Q2,and Q3 internal and normal coordinates are related by a nonlinear transformation. lo Figures 1-4 show perspective views, from several directions, of the envelopes of four representative trajectories where for this system. In terms of the description in section 11, the surfaces were generated with the mesh origin a t the H =. '/2(Pi2+ p22 + P32) + V(Qi,Q2,Q3) (5) potential minimum and relative to a Cartesian axis system and V(Q1,Q2,Q3)is the transformed, modified Sorbie(x,y,z) (6r1,6r2,6a).The solid angle mesh had m = 50, Murrell potential. However, it was felt to be of most to give a total of about 3100 mesh points. Each trajectory interest to study the shapes of the trajectories in the was run for 15 ps spanned by 6 X 104 time steps, with path interpolation a t about 10 points per time step. A mesh of cross wires (rotated in 4 by 4~45' from the viewing (13) R. T. Lawton and M. S. Child, Mol. Phys., 37,1799 (1979);40,773 (1980);44,709 (1981);M. S. Child and R. T. Lawton, Chem. Phys. Lett., direction) has been wrapped around each surface to give 87, 217 (1982). an impression of depth and solidity. (14) J. T. Muckerman, D. W. Noid, and M. S. Child, J. Chem. Phys., In each picture, cross wires have also been drawn 78, 3981 (1983). showing the outline of the three-dimensional potential (15) K. S. Sorbie and J. N. Murrell, Mol. Phys., 29, 1387 (1975).
2742
The Journal of Physical Ch8miStv, Vol. 87, No. 15, 1983
(a>
from
e=
\
56'
q
,
d = 63'
'
'
Ashton and Muckerman
(b)
from
e=
(d)
frcm
8=124'
56'
, d=153'
M
( c > f r o m 8=124O M
, 4=243'
,
$=333'
/
Figure 4. Perspective views of a stochastic trajectory at a total energy of 30 310 cm-'. See Figure 8 for scaling information.
energy contour at the total energy of the trajectory. The surfaces representing these contours were obtained by a . simple modification of the method for the trajectories. Since the molecule has an inversion symmetry in the 6a coordinate, the Sorbie-Murre11 potential has a reflection symmetry about the linear H-O-H geometry, i.e., about the p l h e 6a = 0.4193~.The potential barrier for inversion is at 13017 cm-l above the minimum, and below this energy the potential contours consist of two disconnected wells. This is the case in Figure 1, where only one of the potential wells is shown. Figures 2 and 3 show trajectories at total energies above the inversion barrier, which do not, however, cross the symmetry plane. In each case, the potential contour has been truncated at the linear H-O-H geometry for clarity. The trajectory in Figure 4 exhibits many inversions of the molecule, and it is desirable to represent the motion throughout the whole coordinate space. This was achieved by applying the method of section I1 separately in each half of the potential, generating two abutting sets of surfaces. A similar technique was used to represent the full potential contour. Physically, Figures 1 and 2 show respectively a quasi-
periodic normal mode trajectory and a quasiperiodic local mode (symmetry-breaking) trajectory. The motion is confined to a relatively small volume of the accessible coordinate space and the caustics are well-defined, as evidenced by the relative smoothness of the surfaces. Figures 5 and 6 show cross-sectional views of these trajectories after 60 ps. (After this length of time, the surfaces of minimum distance from the origin have collapsed to a negligible volume, so that the envelopes of the trajectories are represented entirely by the surfaces of maximum distance.) The cross sections in the 6a = 0 plane (Figures 5a and 6a) closely resemble the shapes of corresponding types of trajectory found in earlier two-dimensional studies of this system which omitted (froze out) the bending motion.13J4 Noid, Koszykowski, and Marcus (NKM) have described a practical method for the quantization of three-dimensional quasiperiodic trajectories'l which involves evaluation of the areas of Poincar6 surfaces of section. The condition for an eigentrajectory in a nondegenerate system is 2a(ni + y2)h = Ci$P. dQ (6) where the integrand refers to a specific branch of the
Classical Trajectories in Multidimensional Molecular Systems
0.2
(01 Plane 6a = 0
+
The Journal of Physical Chemistty, Vol. 87, No. 15, 1983 2743
( a ) Plane 6a = O 0.5
0.1 n
h
5
5
N L
N L Y)
IC,
0.0
0.0
-0.1
-0.1
0.0
0.1
6rl ( R )
0.2
e
(b) Plane 6 r , = 6 r 2
0.5
0.1
n
F
c
(c
0
r 2
(b) Plane 6 r q = 0
v)
-3 v)
0.0
d
0.0
d
Lo
rg
-0.1
-0.1
0.0
0.1
0.0
6r, ( A )
0.5 6r2 ( A )
Figure 5. Cross-section views of the normal-mode trajectory shown in Figure 1, together with the total energy potential contour. The interior of the trajectory envelope is shaded. The trajectory was run for 60
Flgure 6. Cross-section views of the locabmode trajectory shown in Figure 2. The trajectory was run for 60 ps.
PS.
curved path along the limiting simple periodic trajectory. Figures 4 and 8 show a stochastic (ergodic) trajectory. The motion defines no caustics, so that the envelope of the trajectory is irregular and ill-defined. If the trajectory were run for a longer time, the envelope would eventually completely fill the potential contour. The NKM quantization method cannot be applied to such trajectories.
semiclassical wave function, the n; (i = 1, 2, 3) are a set of integer quantum numbers, and the Ciare a set of three lines in Q coordinate space defining topologically independent paths through the trajectory. “Topological independence” refers to the manner in which the lines intersect the caustic surfaces. In the present internal coordinate system, Figure 5 shows that suitable paths Cifor this normal mode trajectory would be 6r1 = 6r2,6a = 0; 6rl = -6rz, 6a = 0; and 6rl = 6r2 = 0. For the local mode trajectory, Figure 6 indicates that suitable paths would be 6rl = 0, 6a = -0.05 T ; 6r2 = 0 , 6a = - 0 . 0 5 ~and ; 6rl = 6r, = 0. These lines would in practice have to be transformed back to the normal coordinate space Q in which the kinetic energy of the system is diagonal; alternatively, the proceding analysis could simply be repeated in Q space to choose lines in that space directly. Figures 3 and 7 depict a quasiperiodic normal-mode trajectory whose envelope closely surrounds the limiting pure asymmetric stretch simple periodic traje~tory.’~ For this trajectory, suitable paths for semiclassical quantization would be 6r, = 6rz = 0.16 A; 6rl = 6r2,6a = - 0 . 0 4 ~ ;and a
IV. Conclusion We have described a method for visualizing and examining the envelopes of multidimensional classical trajectories, and have illustrated its utility by application to the coordinate space of a specific three-dimensional system. The method could readily be extended to four or more dimensions by defming an angular mesh of the appropriate dimensionality. More specifically, it could easily be modified to enable visualization of spatially localized quasibound t r a j e c t ~ r i e s , or ~ - ~of the quasiperiodic orbit dividing surfaces which define the adiabatic dynamic barriers in multidimensional reactive ~ c a t t e r i n g . In ~ a three-dimensional (coplanar) reactive system, for example, these dividing surfaces would be two-dimensional curved
Ashton and Muckerman
The Journal of Physical Chemistry, Vol. 87, No. 15, 1983
2744
0.4
-
5
0.2
N
I
(01 Plane 6 a = -0.04
(0)
TI
Plone 6 a = 0
1
L
Lo
0.0
0.2
1
I t
O.O
; -0.2
I
I
I
0.2
0.0
t
0‘ 4
0.5
0.0
6r, ( A )
6r, ( A )
1
(b) Plane 6 r , = 6 r 2
1.0
4
0.5
--
0.0
--
h
t
t i
0
2* c
Y
8
e
-0.2
0.0 6rl ( A )
0.2
-0.5
I
t
0.0 6r, ( A )
0.5
Figure 7. Cross-section views of the normal-mode trajectory shown in Figure 3. The trajectory was run for 60 ps.
Flgure 8. Cross-section views of the stochastic trajectory shown in Figure 4. The trajectory was run for 60 ps.
sheets, typically located in the entrance and exit channels for the reaction. One limitation of the method is that it cannot easily be adapted to the visualization of arbitrary internal structure within a trajectory envelope. Such structure is important, for example, in the semiclassical quantization of quasiperiodic trajectories exhibiting resonance between two or more vibrational modes. In such cases, the internal structure of the caustics must be examined to determine suitable paths for quantization. For these special cases, an alternative method of visualizing the trajectories, which was developed by Noid et al. simultaneously with the present work,16would probably be preferable. The method
of Noid et al. involves “binning” the trajectory path in a multidimensional mesh of cubic cells. Where internal structure is absent or unimportant, however, the present method and that of Noid et al. are probably about equally useful.
Acknowledgment. This research was carried out at Brookhaven National Laboratory under contract with the U.S. Department of Energy and its Office of Basic Energy Sciences. (16)R. A. Marcus, D.W. Noid, and D. M. Wardlaw, J.Phys. Chem., this issue.
