Classical stretching dynamics in methylene - ACS Publications

The Journal of Physical Chemistry, Vol. 83, No. 8, 1979

Classical Stretching Dynamics in Methylene

Sciences, U S . Department of Energy. For 1978-1979 W.H.M. is Miller Professor of the Miller Institute of Basic Research in Science.

983

(5) See, for example, L. I. Schiff, “Quantum Mechanics”, 3rd ed,

McGraw-Hill, New York, 1968, pp 275-279. (6) For a discussion of complex-valued trajectories, see W. H. Miller, Adv. Chem. Phys., 25, 69 (1974). (7) See, for example, (a) K. W. Ford, D. L. Hill, M. Wakano, and J. A. Wheeler, Ann. fbys. (N. Y.) 7, 239 (1959); (b) N, Froman, Ark. Fys., 32, 79 (1966); (c) W. H. Miller, J . Chem. Pbys., 48, 1651 (1968); (d) J. N. L. Connor, Chem. Phys. Lett., 4, 419 (1969). (8) W. Eastes and R. A. Marcus, J . Chem. Pbys., 61, 4301 (1974); D. W. Noid and R. A. Marcus, ibid., 62, 2119 (1975). (9) S. Chapman, B. C. Garrett, and W. H. Miller, J . Chem. Pbys., 64, 502 (1976). (10) I. C. Percival, Adv. Cbem. Pbys., 36, 1 (1977).

References and Motes (1) M. C. Gutzwiller, J . Math. Phys., 12, 343 (1971). (2) W. H. Miller, J . Cham. Pbys., 63, 996 (1975). (3) R. P. Feynman arid A. R. Hibbs, “Quantum Mechanics and Path Integrals”, McGraw-Hill, New York, 1965. (4) See also the discussion by M. V. Berry and K. E. Mount, Rep. Prog. Pbys., 35, 315 (1!372).

Classical Stretching Dynamics in Methylene James L. Rookstool and Christopher A. Parr* Depsittment of Chemistry, The University of Texas at Dallas, Richardson, Texas 75080 (Received August 17, 1978) Publication costs assisted by The Robert A. Welch Foundation

The classical molecular dynamics of a chemically accurate model potential for the 3B1methylene radical are investigated. Natural stretch vibration modes are studied individually at a full range of energies. These modes represent the molecule’s own “most separable” degrees of freedom. The combined stretch dynamics are studied near the zero-point energy. A procedure for obtaining unbiased Monte Carlo initial conditions for mixed-stretch dynamics is outlined. The relation of these classical dynamics to semiclassical quantum numbers is discussed. A pair of local bond modes are discovered in the classical solutions.

I. Introduction In his first paper on Monte Carlo trajectory calculations,l Don Bunker noted the failure of highly vibrationally excited anharmonic, and even harmonic, molecular models to yield sinusoidally varying normal coordinates or conserved normal-modle energies. One consequence of this failure was the impropriety of the use of angle variables, conjugate to the normal-mode actions in the properly random-phase trajectory initialization required for a Monte Carlo study. Instead a less rigorous and less efficient rejection technique’ had to be employed. As the complexity of his models increased beyond triatomic, Bunker and his co-workers developed four more phase sampling techniques,2 each more ingenious than the last, but with his keenly developed sense of the random, he was satisfied with none of them.28 Even for vibrations near zero-point energies (ZPE), Bunker and Chapman3 found the normal-mode description of anharmonic polyatomic (CH,) dynamics to be inexact, though not so flawed that intermode energy transfer rates competed seriously with bimolecular collision interaction times (0.1-0.2 ps). For energies as low as 0,04 ZPE, Chapman and Bunker3 reported normal-mode failure in the form of systematic energy exchange between methyl radical’s degenerate asymmetric bend normal modes. These low energy failures have been found in rigorously harmonic force models as well.4 Again the implication for classical polyatomic dynamics is that neither reactant initialization nor product internal mode energy analysis may be trusted to normal-mode procedures even at low vibration energies where quadratic terms dominate the molecular potential. Bunker compared his near-dissociatively energized triatomic dynamics’ to those investigated by Fermi, Pasta, and Ulam.5 The problem addressed by Fermi et al. was the planar dynamics of a chain with anharmonically bonded links. Instead of the anticipated random nor0022-3654/79/2083-0963$0 1.OO/O

mal-mode energy scrambling, the chain was observed to cycle through a particular sequence of configurations for the entire observation period. The dynamics were locked into a repetitive pattern which excluded an energetically available class of incredibly large measure; thus the phase space of the chain appeared to be metrically decomposable, i.e., significant portions of phase space were dynamically inaccessible to given trajectories. This behavior was in direct violation of the spirit of the ergodic hypothesis, hence vaguely suspect. Bunker noted no such behavior in his dissociative, anharmonic molecular models. “Unless one believes that at some value of (total energy) lower than that used (in his unimolecular decomposition studies) there suddenly appear uninteresting anharmonic molecules”, Bunker concluded,’ “in the molecular problem all points can reach an arbitrarily small volume of phase space and (the Fermi, Pasta, and Ulam) phenomenon is absent”. In his study, “interesting” molecules dissociated, and while a large majority of sufficiently energized harmonic molecules were often “uninteresting” on the time scale of his trajectories (0.8-1.6 ps), all of the anharmonic molecules were “interesting”. Bunker’s study’ and that of Fermi et al.5 clearly lay on opposite sides of what has become known as the “ergodic limit”, a critical energy below which the majority of trajectories are dynamically constrained such that their motion is “regular”, “quasi-periodic”, or metrically decomposable, i.e., Fermi-Pasta-Ulam-like. Above this energy, critical for the suddeness of the ergodic onset, the dominant dynamic behavior is consistent with metric indecomposability; virtually all trajectories can reach an arbitrarily small volume of phase space. This is the character of Bunker’s dissociating anharmonic mo1ecules.l The existence of the ergodic limit for a large class of Hamiltonian problems, including molecular vibration, has been established both theoretically6 and numerically.’ No

...

0 1979 American Chemical Society

J. L. Rookstool and C. A. Parr


964

TABLE I: Properties of t h e ODIM 'B, CH, C,

,C ,

re = 1.082 a f,, = 5.9897 u I = 3361 rn = 1.087 A

B e = 127.7" frr = 0.0768 u 2 = 940 VT*= 9.3 kcal mol-'

D,(HC-H) = 107.5 = 0.0945 u 3 = 3242 cm-' :. ZPEZ 10.8 kcal mol-'

fro

longer suspect, even to statistical mechanics! it encourages the search for nonpathological angle variables for systems below the ergodic limit. It guarantees caustics: boundaries for the turning points of vibrational coordinates, and taxonomic tools for discriminating between isoenergetic systems with different vibration mode energy distributions. For quasi-classical Monte Carlo trajectory studies, not only must vibrational phase (Hamilton-Jacobi) angles be random, but the associated classical actions (hence vibrational energies) must be close to the quantal ones. Indeed several investigatorsg are currently exploring the correlations between classical and quantal actions in multidimensional, semiclassical eigenvalue studies. These studies are predominantly upon model systems of two mathematical degrees of freedom. We report here on our continuing investigation of the classical vibrations of a realistic model of methylene, CH2. Our, as yet unrealized, goal is a proper quasi-classical description of this nonlinear, highly anharmonic radical which we wish to use as a reactant and product in Monte Carlo trajectory studies. While pursuant of a full three-mathematical dimensional method (3N - 6 = 3, asymmetric, bend, and symmetric), we content ourselves here with one- and two-dimensional subsets of the three-dimensional dynamics. In section 11, we review our CH2 potential briefly. In section 111, we discuss some one-dimensional motion, the natural stretching modes, of our molecular model. Finally, we describe the dynamics of the combined stretch vibrations at an internal energy near the model's zero-point energy. 11. 3B1CH, Potential-Energy Hypersurface The potential energies and the associated interatomic forces required for our classical trajectories for ground state (3B1)methylene were obtained directly from an optimized diatomics-in-molecules ODIM calculation.1° The basis set used was the minimum required for full hybridization: hydrogen (%) and carbon (3P,lD, 'S, 5S, and 3D). The optimization parameters (CH VB state mixing function parameters) were set to give the experimental thermochemistry and ODIM potentials resembling those from a careful CI calculation.'l The resulting CH2characteristics are shown in Table I. It is clear that even the zero-point energy (ZPE) is sufficient to surmount the small barrier V,f to the collinear saddle. Although that would require over 85% of the ZPE to collect in the v2 bend mode, a vanishingly unlikely event, it points up the severe anharmonicity of the model's bend potential. (See Figure 2 in ref 10.) The methylene bend mode is thus a serious double minimum problem from the ground state up. Even in the absence of stretch action, the "inversion" is classically allowed for normal-mode quanta u2 I 3. To avoid the problems associated with inversion for this report, we systematically set bond action to zero via techniques described in subsequent sections. Note that this does not imply fixed-angle dynamics. Even in the normal-mode approximation, the non-zero fd ensures that pure stretches will flex the bond angle. The methylene stretching potential is shown in Figure 1. There the five energy contours are given as j z kcal mol-' where j = 1, ..., 5. Were the potential harmonic, V would ~ some orthogonal x , be proportional to (Ax)2+ C ( A Y )for y axis choice, and the dimensions of the contours would

D,(CH) = 83.6 kcal mol-'

f e e = 0.2345 mdyn

a-'

40.18

Figure 1. Equipotential contours for ODIM 3B1 CH, stretches. The vertical is the symmetric stretch coordinate '12(R1 R,); the horizontal is the asymmetric coordinate '12(R,- R 2 ) . The plus sign (+) marks the equilibrium point. The contours are labeled with their energy (in kcal mol-') relative to equilibrium. The bond angle has been set at 126.657° for this figure.

+

increase linearly with j . It is clear that they do not since the contours are not geometrically similar. The potential is visibly anharmonic beyond V = 9 kcal mol-l. Thus dynamics at ZPE = 10.8 or null-bend ZPE = 9.4 kcal mol-' are expected to be nonnormal. 111. Natural Methylene Stretch Modes Individual normal-mode vibrations yield cyclic, rectilinear atomic displacements. Natural vibration modes might not preserve this rectilinearity but must necessarily be periodic. Not every periodic trajectory represents a natural vibration mode, however, as the energy dependence of anharmonic oscillation periods ensures that multipliperiodic, closed trajectories can be found from the beating of natural vibrations. Hence there may be many more periodic vibration patterns than there are independent vibration modes. We seek the vibration patterns which reduce to the normal modes asymptotically at low energies. All trajectories discussed in this paper were obtained from a modified version of the REACS program developed at Toronto.12 Hamilton's equations of motion were integrated by an Adams-Moulton fourth order predictorcorrector method operating a t time steps around s N 0.01 stretch periods. A t CH2 zero-point energies, trajectories could be accurately evolved for at least 45 stretch periods. Such trajectories required 2 h each on an IBM 370/155, but useful results were obtainable from trajectories only 1 / 5 as long. Methylene's symmetric stretch mode is the most easily sought, because the symmetry of the motion automatically reduces the dimensionality of the phase space to searched. Trial initial conditions consist of equivalently stretched CH bonds with bond angle varied trivially about the equilibrium value. When, for a given stretch, the initial angle is found for which the subsequent motion is closed and singly periodic, a natural stretch vibration has been isolated. The total classical action for one such vibration period is thus all symmetric stretch action, and semiclassical symmetric stretch quantum numbers13are readily obtained from these trajectories. Table I1 shows the critical parameters of several such symmetric stretch trajectories. Thus the first excited state is elevated only 3158 cm-l above the ground state, a



965

TABLE 11: Methylene Symmetric Stretch --

UJWKB

-0.5 0.00034 0.7266 0.9996 1.9951 3.0009 4.0004 5.0003 6.0110 8.0052 14.0024

E, kcal mol-’ 0.0 4.609 11.193 13.637 22.398 31.015 39.349 47.456 55.423 70.39 110.5

v,

cm-’ 3242 3200 3141 3119 3038 2918 2876 2795 2716 2592 213,

R, ‘4

e, deg

1.0819 1.0144 0.9806 0.9712 0.9443 0.9237 0.9071 0.8930 0.8806 0.86 0.82

1.0819 1.1611 1.2121 1.2280 1.2791 1.3248 1.3671 1.4074 1.4471 1.52

127.70 127.91 127.98 127.99 128.00 127.96 127.90 127.83 127.74 127.54 126.9

127.70 127.41 127.20 127.13 126.91 126.69 126.48 126.28 126.09 125.74 124.8

discrepancy of 84 cm-l from the normal-mode approximation for this model. From the values of the bond angle at the inner and outer stretch turning points (6, and I%), one sees that rather little angle variation is involved even a t u3 = 14 where the energy is 3 kcal mol-l in excess of the dissociation energy. No natural bend vibration is implied by this variation sincle there is no component of the motion oscillating at or near the bend frequency. Symmetric stretches like these could be found at energies virtually all the way to atomization (191.1 kcal mol-l). They would not be meaningful a t energies above the ergodic limit, which lies well below dissociation, for they would not be representative of the motion the system underwent if perturbed. If one is trying to mimic a quantized molecule, as in quasi-classical studies, the zero-point vibrations of the other two modes would be far more perturbing than necessary to send highly energized “symmetric stretch” oscillating into every energetically available corner of phase space.14 A t energies below the ergodic limit, however, these vibrations might be the “coordinates” of choice for methods which assume vibrational separability. It might be noted that once one is normalized to the fine points of leading one’s target (bend periods are dependent on stretch energies) in this suppression of angular action, a trajectory, closed to seven decimal digits, can be discovered in only seven trial vibration cycles. Thus the search is economical. The natural asymmetric stretch mode has been discussed in ref 15. In this case, the absence of symmetry renders the search for closed, asymmetric trajectories about twice as difficult as for the symmetric ones. The first excited asymmetric state lies 3264 cm-l above its ground state, a 97-cm-’ discrepancy when compared to the normal-mode prediction. Figure 2 shows the natural asymmetric stretch at an energy of 11.193 kcal mol-l, well in excess of the asymnietric zero-point energy of 4.72 kcal mol-l. (This high energy was chosen to entirely encompass total-zero-point dynamics including reasonable perturbations.) Thus the semiclassical quantum number for this trajectory is 0.684. Although bond angle has been projected out in the representation of Figure 2, it varies from 126.23’ at the midpoint to 126.25’ at the end points of the oscillation. The fixed angle chosen for the potential-energy contour, however, is 126.657’. That angle is appropriate to combined symmetric-asymmetric stretch at this energy, which is described in the next section. As can be seen from Figure 2, this discrepancy is irrelevant on the scale of the diagram. Asymmetric stretch is observed to pass outside the equilibrium position, and this curvilinearity becomes very exaggerated at higher energies.15 Contrasting this behavior with the rectilinearity of normal-mode vibration, one sees that the expansion of oscillations in normal modes is an expansion about the wrong center. A similar expansion in the natural vibrations might be expected to be more rapidly convergent. Again, though these natural asym-

1.75

\

Figure 2. A natural asymmetric stretch mode in methylene. The circumscribing curve is a potential curve for 11.193 kcal mol-’. The plus sign (+) represents the stretch equilibrium position. The bond angle for the potential contour is fixed at 126,657’ (see text). The horizontal curve is a closed, periodic asymmetric stretch trajectory at JWKB quantum number vaSym= 0.684 and vasym = 3251 cm-’.

Flgure 3. Perturbed asymmetric stretch. Total energy and potential contour as in Figure 2. Ratio of asymmetric to symmetric energy is

-10.

metric vibrations can be found at very high energies,16they cannot be meaningful guides to the motion at energies in excess of the ergodic limit. Natural symmetric stretch at 11.193 kcal molw1,if plotted on Figure 2, would be a perfectly dull vertical line passing through the equilibrium cross. Its vibration parameters are listed in Table I1 for quantum number 0.7266. That vertical axis and the curvilinear horizontal asymmetric stretch would seem to form a reasonable “vibrational coordinate system” in which to analyze the behavior of general stretch motion at near-zero-point energies. We explore that possibility below.

966

The Journal of Physical Chemistry, Vol. 83,

No. 8, 1979


-

Figure 4. Mixed symmeric/asymmetric stretch. Energy and contour as in Figure 2. Ratio of asymmetric to symmetric energy is 1. Initialization at upper right.

Flgure 5. Same as Figure 4. Initialization at lower right.

IV. Mixed Stretching Vibrations Figure 3 represents a portion of a trajectory wherein the asymmetric vibration of Figure 2 is perturbed by symmetric motion. The total energy is unchanged; indeed all vibrations discussed in this section are at a 11.193 kcal mol-l total energy. The energy partitioning in Figure 3 is about 1 O : l asymmetric:symmetric. For as long as the equations of motion for the trajectory in Figure 3 can be integrated accurately (about 45 oscillations), the vibration remains confined to the region in coordinate space it is shown as occupying in the figure. All turning points lie on the boundaries (caustic curves) of that region. This is characteristic of non-ergodic, quasi-periodic motion. Superposition of Figures 2 and 3 would show the asymmetric stretch of the former almost perfectly bisecting the region of the latter. The natural asymmetric mode is thus stable to perturbations and representative, in many senses, of those perturbed vibrations. Equilibrium lies near the lower caustic; hence normal-mode asymmetric and symmetric stretch, which yield rectangular caustics, would not be good separable coordinates for the motion in Figure 3. However, the natural vibrations of section I11 do not serve this purpose either; neither the top nor bottom caustic have the curvature of the pure asymmetric stretch of Figure 2. Likewise, the side caustics are not vertical. Unless one devised a transformation barrel-distorted on one axis and pin-cushion-distorted along the other, stretch vibration would not be separable in natural vibrations. Thus, although they capture much of the flavor of mixed vibrations, the natural modes are not proper substitutes for the good angle ~ariab1es.l~

Figure 6. Same as Figure 4. Initialization at bottom.

Figure 7. Superposition of curves in Figures 4-6 and their asymmetric images. All six trajectories conform to the same caustic.

Figures 4-6 show portions of three different trajectories wherein asymmetric and symmetric stretch roughly share the available energy. The trajectories in Figures 4 and 5 were initiated with zero kinetic energy from coordinates at two different points on the energy contour. The one in Figure 6 was initiated from a compressed symmetric configuration (bottom center) with only asymmetric kinetic energy. As dissimilar as these trajectory segments might appear, they all display precessional character. They all represent the same distribution of action between symmetric and asymmetric stretch. They are all bound by the same caustic! This can be seen from Figure 7 where the trajectories and their asymmetric mirror images are superimposed. Integration of the equations of motion to the integrator failure point for any of these trajectories continues to preserve the vibrational boundaries shown in Figure 7. Coincidence of these trajectories is even more interesting when it is recalled that figures such as Figure 7 are projections onto the stretch coordinate plane from the full 3N - 6 = 3 dimensional motion. Just as the natural stretch modes flex the bond angle, so do their combinations. Thus the bond angle varies along all six trajectories of Figure 7 , rendering the caustic a surface in a threedimensional space. The trajectories are all embedded in the same surface. The angular depth of the caustic surface can be seen by a view in which the asymmetric coordinate is projected out instead of the bend coordinate. Figure 8 is the object in Figure 7 viewed from its right side, i.e., from the x axis instead of the z axis. Instead of a perspective view of the trajectories themselves, Figure 8 indicates the points at which each trajectory crosses one of eight constant



128O

/ - 7 - 1 - - - T - - - 7

126"

j L

+

I .o

I I

1

1.2

1.1

Flgure 8. Caustic surface of Figure 7 viewed from the asymmetric coordinate [ a = l/z(R1 - R2)]axis. Abscissa is symmetric coordinate [ s = '/hR1 RZ)], and ordinate is bond angle. The plus sign (+) marks equilibrium. Points represent trajectories crossing constant asymmetric coordinate cuts at a = 0.00,0.02,0.04,0.06,0,08,0.09, 0.10, and

+

0.11 A.

I 1

10

I

I

''

' 126'

I

(R,+R,)/Z

(A)

Three-dimensional perspective of the a axis cuts of Figure 8. Front plane is a = 0.1 1 A. Rear plane is a = 0.00 A, Le., the plane through the observer and the symmetry axis in Figure 7. Tic marks on a axis (top left edge of box) show the position of the cut planes. Star in rear plane marks the equilibrium configuration. Figure 9.

asymmetric coordinate planes. The values of the coordinate, defined as a -= ll2(R1- R2),are a = 0.00,0.02,0.04, 0.06, 0.08, 0.09,0.10, and 0.11 A. As can be suspected in Figure 8, and confirmed in the three-dimensional perspective plot Figure 9, the loci of trajectory crossings in a given cut generally form one or more closed curves. Thus the caustic surface is double valued in angle; it forms a surface of at least two intersecting sheets. Figure 8 shows the maximum angular separation of the sheets to be 0.4". On average, the caustic sheets lie about 1' below the equilibrium angle, with extended configurations at smaller bond angles than compressed ones, in keeping with Table 11. Again, though the bond angle varies, no component of any trajectory oscillates with frequencies near the bend frequency or any of its overtones. Bending vibration has thus been rigorously excluded from this problem. The perspective view of the caustic surface, shown in Figure 9, offers further insight into the behavior of near zero-point methylene stretch. The figure shows a > 0. The caustic is symmetric about a = 0; so one can imagine the sheets extending back into the a < 0 direction (not shown) exactly as they extend forward from it. The caustic surface in the a = 0 plane is unique; both sheets have collapsed to a single curve. As asymmetric coordinate a grows through 0.02, 0.04, 0.06, and 0.08 A, the caustic sheets separate, but each such cut indicates that the sheets are pinched together a t some ,point in the cut plane. These lateral intersections are indeed continuous on the caustic surface and represent a one-dimensional seam in the

967

asymmetric direction (with natural asymmetric mode curvilinearity) just as the curve in the a = 0 plane represents a one-dimensional seam in the symmetric direction. Between a = 0.08 and 0.09 A, this asymmetric seam intersects the right-hand caustic of Figure 7 a t its waspwaisted minimum. For the larger values of a , beginning with 0.09 A shown, the surface section has unconnected lobes since the vertical cut through Figure 7 is being made to the right of the waist on the right-side caustic. The cuts a = 0.10 and 0.11 8, are made so far to the right (Figure 7) that only the top lobe is intersected. Trajectories cross this surface (Figure 9) on a top or bottom sheet until the encounter either an edge (the caustics of Figure 7) or a seam (either the symmetric or asymmetric ones described in the last paragraph). Then, as the sheets in these regions join upper and lower surfaces smoothly, the trajectories embedded in them make the transition to the bottom or top sheet just as smoothly. Thus all trajectories heading into or out of caustic corners (Figure 7) are on bottom (smaller angle) sheets whereas those cutting across corners are on top sheets. The reason bottom sheet trajectories heading into corners come out of them on the bottom as well is that they have suffered two edge transitions (top and side, say, for the upper right corner). The only trajectory (zero measure) which suffers only one edge reflection is the one which hits the corner dead-on, landing on the total energy contour. Its kinetic energy goes to zero and it necessarily reflects directly back upon itself, thus preserving the sheet (bottom) upon which it entered. A trajectory streaking from one corner to its opposite, across the middle of the figure, must necessarily cross both seams, thus leaving it on the bottom sheet as it enters the opposite corner. Between the seam crossings, of course, it's a topside trajectory. Of what conceivable utility is this four-blistered pillow? Under circumstances where the suppression of bend action has little chemical or physical relevance,ls the parametrized dynamics represented by the caustic surface enable one to initiate constant-action, mixed-mode trajectories with properly random oscillator phases. This is most easily accomplished by parametrizing the density and (twovalued) "velocities" of trajectories crossing the a = 0 seam as a function of s along the seam. Trajectories initiated from that distribution are but one-fourth of a (parametrized) stretch period away from any arbitrary phase pair. A collision time delay, random over that fractional period, completes an unbiased initiation scheme. A good action-angle analysis," if available, is to be preferred, as it generalizes with less effort to higher dimensionalities. The caustic pillow of Figure 9 yields some suggestions about methods for extracting semiclassical actions in this multidimensional problem. We mention four and carry through with two of them.lg All the methods involve the Einstein-BriIlouin-Keller action integrals on closed C U ~ V C; ~suchS that ~ the ~ desired ~ ~ quantum ~ ~ number ~ ~ ~ is

where 6i = ' I 2for harmonic oscillators and i runs from 1 through 3N - 6 (nonlinear molecule). (1)The two C curves are the symmetric and asymmetric ~ e a m s ,and ~ ~p, and ~ ~ q are taken in projection along the seams. (2) The C curves are symmetric and asymmetric phase space curves generated by trajectories crossing Poincar6 surfaces of sections.ga Those "surfaces" are defined by the seams.

~

~

968



(3) The C curves are the (projected) caustics, and p and q are symmetric or asymmetric projections along the

caustics.20 (4) As in (3), but p and q are taken to lie along the caustics themselves. While (1)and especially (2) seem the most promising, (4) and especially (3) seem to be the most pedagogic. We opt to explore the latter two. The caustics to which we turn our attention are the top and bottom of Figure 7. The difficulty of defining these three-dimensional caustic curves from discrete tangency conditions (from trajectory turning points) cannot be overemphasized. It involves a heavily iterative procedure with diminishing envelopes and several trajectory reruns when older versions fail to interpolate one another to the necessary accuracy (at least lo4 A). It is the fact that trajectories are looping in three dimensions a t the caustics which makes the seams [procedures (1)or (2)] with their near-rectilinear crossings seem desirable. Once the caustics have been defined, local decompositions of kinetic energies at the turning points must be made (see Appendix). The action may then be integrated by any method tolerant of the nonuniform points available from the trajectories. For procedure (3), the interesting test is the action accruing in the asymmetric (A: axis) direction along either the top or bottom caustics. Since the asymmetric action ought notz0depend upon the choice of C, we are surprised that our conservatively error-estimatedz1values for (u + &9ym along the top and bottom are 0.748 f 0.005 and 0.774 h 0.005, respectively. Perhaps, for highly anharmonic caustics, projection along normal-mode axes is inappropriate. For procedure (4), the action along the top and bottom caustics themselves (rather than Cartesian projections) is determined to be 0.813 f 0.005 and 0.808 f 0.005. Thus the discrepancy between asymmetric actions along top and bottom caustics lies comfortably within the generous error bars. From this procedure the semiclassical asymmetric quantum number appropriate to Figure 7 is judged to be about 0.31.

V. Methylene Bond Modes In the last section, we confirmed and explored the stability of methylene's asymmetric stretch with respect to symmetric perturbation at zero-point energies. We observed caustics to grow steadily away from the natural asymmtric mode as the fraction of symmetric energy increased from zero to near 50%. This asymmetric stability in methylene is contrasted by its symmetric instability. While, as indicated in section 111, it is possible to find pure natural symmetric stretch trajectories at all energies from zero to atomization, this mode does not characterize the motion resulting from even trivial asymmetric perturbations at the near-zero-point energy of section IV (11.192 kcal mol-l). The behavior will be demonstrated in Figure 10 by a huge 0.3% perturbation, so that the motion will be visibly apparent. The trajectory therein was initiated with zero kinetic energy, on the potential contour with slight positive asymmetry. The motion begins as near-vertical (symmetric) oscillations which diverge exponentially from that axis in a libration which reverses on the tilted portion of the caustics and repeats with a cycle time which becomes arbitrarily long as the perturbation is diminished. Negative perturbations show an identical but mirror-imaged caustic. The motion for a perturbation many orders of magnitude smaller would display virtually the same caustic as shown in Figure 10. The only difference would lie in the high density of vertical strokes which would render virtually black the narrow

Figure 10. Perturbed symmetric stretch. Contour and energy as in Figure 2 . Ratio of asymmetric to symmetric energy is -0.003.

Flgure 11. Perturbed bond mode. Contour and energy as it7 Figure 2 . The plus sign (+) marks equilibrium.

vertical strip between the symmetry axis (note the tic marks in Figure 10) and the currently near-vertical caustic segments. The trajectory in Figure 10 represents an example of one which will cyclically diverge and converge exponentially upon a perturbation of itself.22 While, except in a time-averaged sense, symmetric stretch does not accurately reflect its perturbed motion, there is a stable mode of motion at this energy which does capture the essence of this caustic. As the asymmetric perturbation grows, the caustic shown in Figure 10 does not expand; rather it collapses, eventually to a closed, periodic trajectory, which is, itself, stable to stretch perturbations. A small perturbation on this mode is shown in Figure 11 to indicate both its stability and its orientation. Just as unperturbed symmetric stretch implies equal amplitude, in-phase bond motion, and unperturbed asymmetric stretch implies equal amplitude, out-of-phase bond motions, so too the mode of Figure 11yields in-phase bond motion but with unequal and periodically varying bond amplitudes. Due to the relatively heavy mass of the central atom, the bond vibration axes are very nearly at 45 and 135' in Figure 11. From an analysis of the end points of the mode depicted in that figure, one finds that the bond amplitude ratio AR2/AR1 varies from about 0.26 at the compressed to 0.33 at the expanded configurations. That implies a bond potential-energy ratio variance from about 0.07-0.11 over the same range. At the vibration midpoint, the trajectory again yields dR2/dR1 = 0.26 for a bond kinetic-energy ratio of about 0.07. Thus the mode in Figure 11 is oscillating with, over 90% of its energy in R1. It could thus justifiably be called a local bond mode. It is a natural consequence of the classical motion at this


energy. However, the action ratio is such for this vibration that 0 < u, < 1 while -0.5 < u, < 0; hence it should have no quantum correspondent. If this stable bond mode persists at higher energies, an energy can be found where its caustics support quantizied actions. At that point, the nature of the exact quantum solution would be very interesting. To conclude our discussion of the low-energy stretches of methylene, we must provide the transition behavior between caustics like Figure 7 and those like Figure 10. The transition occurs as the midpoints of the top and bottom caustic of Figure 7 arch toward their corresponding points on the contour as symmetric energy fraction increases. Simultaneously the caustic corners move farther apart until, for example, the lower left corner of Figure 7 is nearly coincident with the lower left corner of Figure 10, and likewise for their right corners. At that energy distribution, Figure '1's top and bottom caustics are tangent to the contour, forming there caustic corners like the near vertical ones of Figure 10. The expanded Figure 7-like caustics then splits into a symmetric pair of Figure 10-like caustics, and the transition has occurred. An elegant description of a similar phenomenon in a Csusystem is given in ref 9a. Acknowledgment. We are grateful to Professor J. Moser for helpful discussions. We acknowledge the generosity of The University of Texas at Dallas' School of Management and Administrative Sciences and the Office of the Vice President for 13usiness Affairs for their donation of computer time for this project. This work was supported in part by The Robert A. Welch Foundation and the Organized Research Fund of The University of Texas System. We thank Schatz and Moser for a preprint of ref 17.

Appendix. Local Separation of Action Components A dynamical system evolving along its trajectory accumulates action A at a rate equal to twice the kinetic energy: dA/dt = 2% dA = 2T d t = 2T(dt/dql)dql

If there exists a local orthogonal coordinate system then

T=

1/22Cr1Qz2 1

which allows the time integration to be changed to coordinate integration dA = 2 X j/2CyLQ12(dt/dql)dq]= I

where (dq;)/(dql) is a measure of direction of the trajectory near the caustic rather then (dq,/dq,). The mass weighting factors y i can be evaluated by (1)fixing coordinates qz (i # j ) , (2) transforming an elemental displacement in ql into (center of mass) Cartesian coordinate displacements, (3) obtaining the associated Cartesian kinetic energy, and (4) dividing by Y,~. Alternately, the Cartesian coordinates, expressed as functions of the trajectory's own (internal) integration coordinates, can be differentiated with time; the resulting kinetic energy expression is rearranged into quadratic form in the internal coordinates. Our internal coordinates are s = 1/2(R1+ R2),a = llz(R1 - Rz),and 0. At fixed 0

The Journal of Physical Chemlstty, Vol. 83, No. 8, 1979

T' = 2mHh2[1- B cos2 (a - 0/2)]

969

+

2mHs2[1 - B sin2(a - 0/2)]

+

where B = 1 - [ m ~ ( 3 m + ~mc)/M] and M = mc 2mH, the total mass. Although the expression for To is more complicated, it is completely dominated, for our case, by the term ys02where yo = 2 m ~ / ( R 1 -+~ R2-2) Thus "asymmetric" kinetic energy becomes T-~ 0 8 ~ T, = s2(1- B sin2 0 ) 1+ h2(1 - B cos2 8) (A similar expression obtains for T,.) This form allows the kinetic energy along a curve C at point P to be obtained from the potential at P, the total energy and the squared trajectory slope (ds)2/(da)2. Care must be taken to interpolate slopes for an arbitrary C along the same branch of the trajectoryaga

References and Notes D. L. Bunker, J. Chem. Phys., 37, 393 (1962). (a) D. L. Bunker and W. L. Hase, J . Chem. Phys., 59, 4621 (1973); (b) H. H. Harris and D. L. Bunker, Chem. Phys. Lett., 11,433 (1971). S. Chapman and D. L. Bunker, J . Chem. Phys., 62, 2890 (1975). C. A. Parr, A. Kuppermann, and R. N. Porter, J . Chem. Phys., 66, 2914 (1977). E. Fermi, J. Pasta, and S. Ulam, "Studies on Non-Linear Problems':. I, LA-1940, Los Alamos Scientific Laboratory, Nov 2, 1955. This study is also reported in E. Fermi, "Collected Papers", Vol. 11, Chicago University Press, Chicago, 1965, p 978. (a) A. N. Kolmogoroff, see Appendix D of R. Abraham, "Foundations of Mechanics", W. A. Benjamin, New York, 1967; (b) V. I. Arnold and A. Avez, "Ergodic Problems in Classical Mechanics", W. A. Benjamin, New York, 1968; (c) J. Moser, "Stable and Random Motion in Dynamical Systems with Special Emphasis on Celestial Mechanics", Princeton University Press, Princeton, 1973. (d) A readable account of what is known as KAM theory (Kolmogoroff-Amokl-Moser) appears as Appendix A.3 in R. Balescu, "Equilibrium and Non-equilibrium Statistical Mechanics", Wiley, New York, 1975. J. Ford, Adv. Chem. Phys., 24, 155 (1973). I. Prigogine, A. Grecos, and C. George, Celes. Mech., 16, 489 (1977). E.g., (a) D. W. Noid and R. A. Marcus, J. Chem. phys., 67,559 (1977); (b) N. C. Handy, S. M. Colwell, and W. H. Miller, Faraday Discuss. Chem. Soc., 62, 29 (1977); (c) I. C. Percival and N. Pomphrey, Mol. Phys., 35, 649 (1978); (d) K. S. Sorbie and N. C. Handy, ibid., 35, 1319 (1977); (e) and references in the foregoing. C. W. Eaker and C. A. Parr, J. Chem. Phys., 64, 1322 (1976). V. Staemmler, Theor. Chim. Acta, 31, 49 (1973). Some other trajectory calculations using this program are reviewed by J. C. Polanyi and J. L. Schreiber in "Physical Chemistry-An Advanced Treatise", Vol. VIA, "Kinetics of Gas Reactions", H. Eyring, W. Jost, and D. Henderson, Ed., Academic Press, New York, 1974, Chapter 6, p 383. See, for example, A. Messiah, "Quantum Mechanics", Vol. 1, North-Holland Publishing Co., Amsterdam, 1970, pp 239-241. (a) P. Brumer and J. W. Duff, J. Chem. Phys., 65, 3566 (1976); (b) E. Thiele and D. J. Wilson, ibid., 35, 1256 (1961). C. A. Parr and J. L. Rookstool, ACS Symp. Ser., No. 58, 250 (1977). We were still locating them above 90% of the atomization energy when we lost interest. The curvilinearity tends toward 90' as these trajectories cross the atomization plateau. This leads us to suspect they persist at energies in excess of atomization, bouncing back and forth between and two bond valleys. For a derivation of good angle variables, see G. C. Schatz and M. D. Moser, Chem. Phys., 35, 239 (1978). Since the bend frequency is so very much smaller than the stretch frequencies, it is not difficult to imagine stretch and bend dynamics at zero-point energies to be adiabatically decoupled. That cannot, of course, be rigorously true, but the error incurred in that assumption might be acceptably small. The difficulty of calculating three-space curves lying on the caustic surface (and not being trajectories), and the sens actions to the accuracy of these curves, dictate the more modest goal. W. Eastes and R. A. Marcus, J . Chem. Phys., 81, 4301 (1974). Caustic curves refined to stationary actions. Aiternate crossing points deleted to test the trapezoid rule integration. Resultant enor estimates doubled. J. W. Duff and P. Brumer, J . Chem. Phys., 67, 4898 (1977).

Classical stretching dynamics in methylene - ACS Publications

Recommend Documents