Defective modelling of chaotic motions on empirical potential energy surfaces

The progression from regular to chaotic behaviour with increasing vibrational energy is examined for two unimolecular reactions, CH 3 NC = CH 3 CN and NCNC = NCCN. The potential energy surfaces used are, respectively, a piecewise empirical construction and an ab initio numerical surface. It is found that in the former case, the motions never become chaotic in the appropriate energy range, but in the latter, they seem to be approaching that ideal condition. The reasons for this difference are subject to speculation at the present time, but there seems to be a strong impediment to randomisation of energy in one case that is not present in the other. An attempt is made to formulate a semi-quantitative measure of chaotic behaviour in these reactions. Until this problem with synthetic potential energy surfaces can be resolved, these results have important consequences for the numerical modelling of larger polyatomic systems, up to and including such problems as protein folding.


Introduction
A cornerstone of the theory of unimolecular reactions is that in the reactive energy range, the internal motions of the atoms are random, or that all energy states are equally probable, or some other statement implying totally chaotic motion. A second cornerstone is that states depleted by reaction are replaced on a time scale much faster than the decay processes; in the case of a single molecule, or at the low-pressure limit, this randomisation process must be first-order, although it may have a second-order component at higher pressures. 1,2 The mechanism of this process may be envisaged as follows: we start from a ground-state molecule, for which the internal motions are regular and vibrational spectra can be observed, but as the energy rises, due to anharmonicities, a pair of spectral lines will merge, i.e. have the same energy, and the wave function becomes indeterminate. The higher the energy, the more crossings there will be, and the degree of chaos increases. For small (triatomic) molecules, there may be no crossings, or else very few, and so their reactions do not fall within the realm of traditional (RRKM) unimolecular reaction theory; on the other hand, large molecules will have many crossings and readily fulfil the requirement of chaotic motion at reaction energies. Of the two isomerisation reactions chosen for this study, the first is generally regarded as falling into the large-molecule class, whereas the other, containing only four atoms, has to be regarded as a borderline example.
In previous studies, 3 we had observed that for the isomerisation reaction NCNC = NCCN, the calculated lifetimes obeyed microscopic reversibility, i.e. r(E) CN /t CN = r(E) NC /t NC where t is the ensemble average lifetime and r(E) is the density of states at the given energy E; however for the reaction CH 3 NC = CH 3 CN, this was not true, with a discrepancy of about an order of magnitude, but improving slowly with increasing energy. The study presented here attempts to understand this departure from one of the fundamental principles of statistical mechanics.

Computational procedure
The computational methods used have been described extensively in previous publications: 3,4 briefly, to follow each trajectory by using a simple fourth-order Runge-Kutta start-up procedure followed by a fifth-order Adams-Mouton predictorcorrector. 5 For a given value of the energy, trajectories were started from the equilibrium internuclear configuration of either NCNC or NCCN with quasi-random momenta assigned to each of the 12 degrees of freedom (or 18 for the CH 3 NC or CH 3 CN cases); a portable random number generator (ran1. for) 5 was used to create these random initial conditions. A time step-length of 0.1 fs was used throughout, and values of all internuclear distances for the 4-and 6-atom reactions were tabulated at dt = 5 fs intervals.
Each of the (r 2 r e ) N ði;tÞ vectors for trajectory N, i = 1,i max , t = 1,t max , (i max = 14 and t max = 16 384dt in this study -see below), was subjected to a Fourier transform algorithm 5 to extract the characteristic frequencies; however, the individual frequencies themselves are of minor importance because, separately, they contain no angle information, but they do provide an intuitive feel for the connections of specific motions with particular spectral frequencies.
It is a simple matter, for the 4-atom case, to see that a vector containing the sum of all 6 distances at each time step contains all the information, internuclear distances, bond and dihedral angles, and gives the full spectrum.
The 6-atom case is more complicated: there are 15 interatomic distances, but there are only (3n 2 6), i.e. 12 internal degrees of freedom when n = 6. Thus, by using all 15 distances, the system would appear to be over-specified, but not so: it is feasible, but not trivial, to transform these 15 vectors into 12 vectors, each representing one normal mode, but complicated further because there are 4 doubly-degenerate as well as 4 nondegenerate modes. Moreover, the concept of a normal mode only has meaning for low energies where the motions are regular, and it is unclear how instructive such calculations would be when applied to partially or completely chaotic motions, so this avenue was not explored further.
A useful contraction is to impose the constraint that all three H atoms are grouped with one C atom to give a pseudomolecule N-C 1 -C 2 -M where M represents the centre of mass of the three H atoms; as before, N-C 1 -C 2 -M has 6 interparticle distances, and the 4-atom group C 2 H 3 contains another 6, making 12; however, the rotation of the CH 3 group when N-C 1 -C 2 -M is non-linear is not defined, and it needs distances from any one of the H atoms to C 1 and to N to specify the dihedral angles completely, making 14 in all. It was thought that this parallelism between the two systems might be of help in their comparison.
The calculations are subject to two constraints: (a) that the trajectory shall be long enough to contain a large number of vibrational periods, and (b) that the sampling width for the Fourier transform (FFT) algorithm shall be wide enough to give a good resolution. With time intervals of 5 fs, a choice of 256-point sampling width gave reasonable results; this leads to a 16 384-length vector for each variable, and a time length of 81.92 ps. At energies near the higher limit of these calculations, there are often many more than 50 crossings and recrossings of the energy barrier on this time scale, 6 whereas only trajectories with no (or very few) crossings are suitable for this experiment; for example, in one of the two sets of high-energy (7.3 eV) calculations, out of 1000 trajectories not one lives for this length of time without a crossing, and only one at 7.0 eV.
The potential energy surfaces used were: for CH 3 NC = CH 3 CN, the empirical surface constructed by Sumpter and Thompson; 7 and for NCNC = NCCN, an 870-point ''clampednuclei'' MP2/6-31G* ab initio potential surface. 6,8 Fig . 1 shows the key properties of the two potential energy surfaces. Since the energy difference between the minima of the CH 3 NC and CH 3 CN potentials is almost exactly 1 eV, it is convenient to use electron volts in these calculations: 1 eV = 23.06 kcal mol 21 = 96.5 kJ mol 21 = 8066 cm 21 ; in what follows (except for Fig. 5) energies are uniformly measured from the global minimum, i.e. NCCN or CH 3 CN. The zero-point energy (zpe) for either CH 3 NC or CH 3 CN is 1.3 eV, so that the highest energy, 7.3 eV, is 5 eV above the CH 3 NC zpe, i.e. 115 kcal mol 21 , which is far above the activation of 38.2 kcal mol 21 . 9 Likewise, the NCNC minimum is about 1.38 eV above the global minimum and the zero-point energies for both NCNC and NCCN are 0.42 eV; this makes the topmost energy tested of 4 eV approximately 2.2 eV above the NCNC zpe.

Numerical results: CH 3 NC and CH 3 CN
At energies near the zpe, the resolved spectral lines are very sharp. but they grow in width with increasing energy. Fig. 2 shows the degree of variation for two individual Fourier transforms, (a) for R(N-C) and (b) for R(C-C). There is a slight shift in frequency of the single peak in the former case, but no evidence of contamination from other frequencies in the system; in the other, in addition to the broadening, there is an apparent shift of about 400 cm 21 of one peak to higher frequencies. This implies that there is virtually no significant coupling of R(N-C) with any of the other motions in the molecule, even when the total internal energy is 5 eV greater. For the R(C-C) case, there is evidence of anharmonicity, but at higher energy, one would expect the energy levels to be closer together, and not further apart. 3  R(C-H) behaviour, and (b), that of the R(H-H) motions. These spectra show much more evidence of mixing and chaotic behaviour, and of frequency shifts in the expected direction. Fig. 4 compares the total spectra of two trajectories. For the 2.3 eV case, these spectra are very similar for each individual trajectory, but vary in the observed intensities due to the differences in starting momenta in each calculation. However, for the higher energy spectra, these are almost identical for trajectories with no (or only a few) CH 3 NC = CH 3 CN crossings; since each of these trajectories starts from a different set of momenta, this sameness suggests comparable degrees of chaos in two isolated regions of phase space.
A semi-quantitative measure of the degree of chaos?
If the extent of the chaos revealed by a spectrum such as Fig. 4(b) were perfect, the spectrum would be a flat line over the whole range, with the total signal strength residing at zero frequency. This plot shows some of the required characteristics: there is a long low signal between 200 cm 21 and 800 cm 21 , and a longer one between 1200 cm 21 and 2000 cm 21 , which could be characterised as something analogous to ''white light''; also, there is a hint of the expected rise near zero frequency.
However, there are ''blockages'' between these two regions of ''almost-chaos'' -a large one near 1000 cm 21 and a weaker one near 1300 cm 21 , which are frequencies associated with R(C-C) and R(H-H) motions in Fig. 2 and 3. At the present time, the exact cause of these impediments to the onset of chaos throughout the system are yet to be understood.
Examining the complete spectra extracted from non-crossing trajectories for all four of the molecules studied in this paper, it appears that there is a relationship between the total area under the curve and the energy: this is shown for CH 3 NC in Fig. 5(a). Even though there is a large scatter in the data, it is clear that there is something like a linear correlation between the area and the energy; the area units are arbitrary, as each point is simply the sum of the 256 values of the strength given by the Fourier transform algorithm for each individual trajectory. The straight line is simply a guess, drawn by eye, and it does not start at the origin, as it probably should; in fact, at very low energies (e.g. in the classical region below the zero-point energy), the motions are too feeble for the Fourier transform algorithm to extract any frequency peaks, due to insufficient arithmetic precision. The large scatter also arises from a combination of the limitation in the trajectory length and in the sampling width, whence the whole frequency range  from zero to slightly over 3300 cm 21 is covered with only 256 points.
As energy is removed from the regular to the chaotic regime, the area under the regular spectrum will decline, and be transferred to the feature near zero frequency. To test this hypothesis, an arbitrary division was chosen at 250 cm 21 (just under the lowest bending frequency of 263 cm 21 in CH 3 NC) and the area below that was compared with the total. For both CH 3 NC and CH 3 CN trajectories having no or few crossings, values of less than 0.1% were found at the lowest energies and up to around 10% at very high energy. This is shown in Fig. 5(b) which exhibits similar scatter, and for the same reasons.
If it were possible to reduce the scatter of these plots, a more accurate relationship between the energy and the area under the spectrum would ensue, and the percentage of the total area residing below some cut-off frequency would yield an acceptable measure of the chaos in the internal motions for each trajectory. This is not quite as straightforward as it may seem: first, it would be necessary to use a larger sampling width, say 512 or 1024, which would necessitate the construc-tion of, respectively, (14 6 32 768) or (14 6 65 536) matrices of interparticle distances as a function of time. This cannot be done by extending the time span because already (see above) very few trajectories of higher energy live for 81.92 ps without a crossing, so the time intervals must be reduced. That will make the inter-atomic differences at each step smaller, with loss of precision, and may require quadruple (32 decimal place) precision, especially at the lower energies.

Numerical results: NCNC and NCCN
The results for trajectories of these molecules are quite similar to those described above at low energies, as shown in Fig. 6; if we label the atoms in this molecule as N 1 C 1 N 2 C 2 , then panel (a) represents the transform on R(C 1 -C 2 ) and (b) on R(N 1 -C 2 ). For NCCN at 2.0 eV, the average percentage of the total area for the 0-250 cm 21 range in 13 trajectories having zero or one crossing is 10.8%, and for NCNC over 11 low-crossing trajectories, 24.6%. At energies greater than 3.0 eV above the NCCN minimum, there are very few low-crossing trajectories Fig. 4 Complete spectral plots, i.e. Fourier transforms of Sr N ði;jÞ , at energies of (a) zpe and of (b) (zpe + 5 eV); in (a), for this particular trajectory, the doublet near 3000 cm 21 is very prominent. in an ensemble of 100, with percentage for the 0-250 cm 21 range being (for 2 trajectories in each case) 13% and 25%, respectively. An example of one of these spectra is shown in Fig. 7(a): the signals around 1000 and 2000 cm 21 from Fig. 6 still persist, but are hardly perceptible on the scale of this diagram: the feature between 400 and 800 cm 21 is probably the remnants of a bending process. 10 Usually, there are very many crossings at these energies for each trajectory, typically in the range 50-60, 6 and the outcomes are quite different. For example, Fig. 7(b) shows the result for a trajectory with 53 crossings over the 81.92 ps period: here, one finds a truly random result -with in excess of 99% of the signal strength concentrated below 250 cm 21 . Moreover, the feature near 1000 cm 21 is probably associated with the crossing frequency: of the 53 crossings, 15 are spread fairly evenly between 9 and 11 ps, i.e. mean separation of 0.13 ps; two other bursts of crossings have mean lifetimes of 0.17 ps, hence, the mean there-and-back crossing time is about 0.3 ps, corresponding to a frequency near 1000 cm 21 .
Earlier spectra purporting to show ergodicity, e.g. Noid, Koszykowski and Marcus, 11 are remote from the asymptotic behaviour of Fig. 7(b) here.

Discussion
The results of this study show conclusively that the synthetic potential surface for the CH 3 NC = CH 3 CN reaction in rotational state J x = J y = J z = 0 is deficient because there is some obstacle to randomisation. Calculations with non-zero J are very little different because the rotational broadening of the lines is minimal, and the rotations don't add significantly to the probability of eigenvalue crossings. The simplest way to describe it would be to say that if the three R(C-H) spectra for J x = J y = J z = 0, one of which is shown in Fig. 3(a), and the corresponding three for (say) J x = J y = J z = 10 for the same energy and initial conditions, were laid out on the desk, one could not guess which ones went with which J; likewise for the three R(H-H) spectra, one of which is shown in Fig. 3(b).
As we have noted previously, 3,4,12 the calculated rate constant for the isomerisation of CH 3 NC is about a factor of 20 below the observed reaction rate, and a similar discrepancy of a factor of 20 was found by Stimac and Barker 13 during the simulation of the CH 3 -O-ON-O radical thermal dissociation into CH 3 O + NO 2 , using a similarly constructed potential energy surface. In addition, their ensemble decay resembles  that of a small molecule like HNC, 14 instead of a first-order process like that for CH 3 NC. 15 That these deficits are due to insufficient mixing has been demonstrated by adding more coupling between the vibrations: this was done by either using coupling constants published by Sumpter and Thompson,7 or by an artificial coupling procedure in which some momenta were reversed every 0.1 ps 12 mimicking, crudely, collisions at very high pressure, or black-body radiation 3 in a collisionless environment. However, there is no way of knowing, beforehand, how much coupling to add, as the reaction rate can be driven too high, especially just above threshold, as shown in Fig. 6 of ref. 12.
Beyond that, there are other issues that need to be resolved. For example, in the light of Fig. 4(b), why does the decay of an ensemble of 1000 trajectories behave in a perfectly logarithmic fashion 15 when there are two isolated quasi-chaotic regions of phase space? This could only be true if molecules trapped on either side of the blockage decayed with the same intrinsic rate.
This raises the question how much chaos is needed for there to be a first-order decay? Both CH 3 NC and NCNC exhibit first-order ensemble decay, 8,15 but in the first case, the distribution of states is far from chaotic, whereas in the second it seems to be approaching total chaos: thus, first-order decay does not imply a proper chaotic distribution among the reactant states. On the other hand, microscopic reversibility, which fails badly in the first case, but not in the second, seems to be a more reliable indicator.
Finally, does microscopic reversibility guarantee the calculation of the correct rate constant? We know that for the isomerisation of CH 3 NC where there is a discrepancy of about an order of magnitude, 4 the calculation of the rate fails badly; unfortunately, there are no experimental data available for the NCNC isomerisation, as it polymerises readily at 193 K, 16 so we may never know.
Recently, Sutcliffe and Woolley 17 have pointed out that clamped-nuclei potential energy surfaces are inconsistent with quantum mechanics and cannot give a correct representation of the transition state problem; however, the Eyring-Polanyi concept of the transition state as a saddle-point on a potential energy surface has been a powerful visualisation tool, and one would expect the numerical errors to be comparable with those usually assosciated with the Born-Oppenheimer approximation. Thus, it is important to determine how well the classical treatment of motion on a clamped-nuclei potential energy surface can match the experimental behaviour.
The smallest molecules for which extensive experimental isomerisation data are available are methyl isocyanide 18,19 and cyclopropane, 20 6-atom and 9-atom molecules respectively. Calculation of ab initio clamped-nuclei potential energy surfaces for the former reaction, perhaps even the latter, are within the bounds of modern computing power, 21 whence it would be possible to test the predictive accuracy of two methods: first, the classical trajectory method used here; and second, a wave-packet method 22 which is essentially a wave-mechanical version of the same classical method.{ The possibility of an ''exact'' solution, based on Sutcliffe and Woolley's analysis, seems to be further into the future.
Finally, from time to time, Ehrenfest's Theorem and/or the correspondence principle are mentioned in connection with such trajectory calculations. In the quantum sense, this association is meaningless because as the energy rises and the first eigenvalue crossing occurs, the adiabatic principle is violated. But the present calculations show that such crossings are nowhere nearly as frequent as one might expect, and often the ground-state vibrational frequencies persist up to very high energies. To resolve this paradox, it should be possible to construct an arbitrary 3-atom potential energy surface with anharmonicities chosen in such a way that the vibrational frequencies will cross at some energy. At this point, what happens? Do they become one more complex motion, is there an avoided crossing, or do they simply cross and continue unchanged?
If either of the last two, by analogy with the high-energy behaviour of a harmonic oscillator, 24 success may be interpreted as a multi-dimensional extension of the correspondence principle 25 with a collection of weakly dependent anharmonic oscillators whose wave functions are similar to those of a simple harmonic oscillator, but with asymmetric distortions. 26

Conclusion
Classical trajectory calculations for two unimolecular isomerisation reactions, NCNC = NCCN using a rudimentary ab initio potential energy surface, and CH 3 NC = CH 3 CN using an empirical potential energy surface, have been examined. These surfaces are quite primitive, but exhibit correct asymptotic behaviour and approximately correct energy values for the stationary states. No effort was made to calculate rate constants, only to examine the difference in behaviour with respect to microscopic reversibility between the two different basic types of potential energy surface.
It was shown, convincingly, that the primitive ab initio surface allowed the atomic motions to become almost chaotic, as required by unimolecular reaction theory, but that the empirical surface did not. The details of this failure were not revealed by these calculations, but (quoting Bowman et al. 21 ) ''simple functional forms are highly problematic'' and it may be a consequence of the use of switching functions in the empirical surface.
If this difficulty can be resolved, and if classical trajectory calculations can give satisfactory reaction rates for one, or more, of the reactions 17-20,23 noted above, it could reveal how to create semi-empirical potential energy surfaces that would allow the prediction of useful rate constants for complex { A third option might be the cis-trans isomerisation of 1,2-dideuterioethylene where, due to symmetry, generation of the ab initio potential surface is simpler, but the experimental data 23 are less extensive. reactions such as occur in combustion, or in biological systems.