Validating potential energy surfaces for classical trajectory calculations

Two possible methods for probing ﬂ aws in empirical potential energy surfaces are examined, using a simple unimolecular isomerisation as an example. One is a heuristic estimate of the degree of chaos in a vibrating molecule; the other is the energy distribution pro ﬁ les at the instant of reaction. It is found that deliberate increases in the degree of computational chaos do not alter the ensemble average of intrinsic vibrational chaos in these trajectories.


Introduction
Since the advent of Eyring-Polanyi transition state theory, it has seemed almost axiomatic that a classical trajectory calculation for a unimolecular reaction on a carefully constructed potential energy (PE) surface could yield an acceptable reaction rate; this assumption rests primarily on the notion that classical and quantum properties converge at high energy, i.e. the correspondence principle will hold for the vibrational motions of a polyatomic molecule, which seems to be a reasonable assumption. 1 The discussion here will be conned to unimolecular reactions, but these considerations apply equally to bimolecular reactions, and beyond. 2 Hundreds of thermal unimolecular reactions have been reported over the past 60 years; 3 among those, three of the most extensively studied, the isomerisations of methyl isocyanide, 4,5 of cyclopropane, 6 and the thermal dissociation of di-tert-butyl peroxide, 7 present clear challenges for computational modelling.
At the present time, ab initio potential energy surfaces are not available for any of these reactions, but for the simplest case, CH 3 NC # CH 3 CN, one can be expected before long.For the other two target reactions, an ab initio potential energy surface for the isomerisation of cyclopropane to propylene may follow later, but probably not for the thermal dissociation of ditert-butyl peroxide, for which one would need to create an empirical surface.But, as Bowman et al. 8 have said, empirical surfaces are "highly problematic", and the purpose of the work described below is to probe possible methods for exploring deciencies.
In fact, such deciencies arise in the simplest of these three cases, that of classical trajectory calculations for the reaction CH 3 NC # CH 3 CN using the well-known empirical potential energy surface of Sumpter and Thompson. 9There, the asymptotic forms at large and small interatomic distances are realistic, and the vibrational properties near the equilibrium structures are well represented.2][13] As TST is an equilibrium theory, 14 it cannot illuminate dynamical properties of the reacting molecules in ways that trajectory methods do, and so these conicts need to be understood.
An earlier study 12 had conrmed that a similar calculation for the reaction NCNC # NCCN, but using an ab initio potential energy surface, gave results consistent with microscopic reversibility, implying that the CH 3 NC # CH 3 CN potential energy surface must be defective somewhere in an intermediate energy region.It had long been suspected that the cause may be a lack of sufficient vibrational mixing in this energy range, 11 and led to my development and initial testing of a speculative method for quantifying the degree of chaos in a vibrating molecule. 13This involved the derivation of a "mechanical" spectrum as the Fourier transform of the sum of all inter-atomic distances in the vibrating molecule; this type of spectrum should not be confused with a conventional experimental electromagnetic spectrum.Some features of such a spectrum could be identied by comparison with transforms of individual distances, or of groups of distances, and proved to be helpful in the analysis of the present calculations.
Most of the emphasis of this paper is focussed on the CH 3 NC molecule, but this potential energy surface has its origin at the potential minimum of CH 3 CN.As noted previously, 13 the electron-volt (1 eV ¼ 23.06 kcal mol À1 ¼ 96.5 kJ mol À1 ¼ 8066 cm À1 ) is a useful energy unit as it puts the potential minimum of CH 3 NC at 1.0 eV with respect to the global minimum.The zero-point energy (zpe ¼ 9626 cm À1 ¼ 1.19 eV) of CH 3 NC then stands at 2.1980 F 2.2 eV, and the potential barrier at $2.8 eV; the zero-point energy of CH 3 CN (9623 cm À1 ) is also 1.19 eV, but is not used in the present work.
Anchoring the energy measurements to the minimum of the CH 3 CN system can lead to confusion, but I have tried to make this distinction clear all through.Energies in eV are measured from the potential minimum of CH 3 CN unless otherwise stated, with an explanatory comment if needed; wavenumbers will be used mainly for state densities and spectral properties, and values imported from elsewhere are cited in their original units.
Table 1 lists a selection of relevant data for the forward reaction CH 3 NC / CH 3 CN from earlier calculations, 11 to provide some justication for the classical approach.It shows that a quantum state having any of these lifetimes would overlap with thousands of neighbouring states, 15 so that its trajectory may explore the whole of the available phase space within that energy band.These energies start somewhat above reaction threshold (the experimental activation energy 16 is 38.2 kcal mol À1 , i.e. 1.66 eV above zpe, 3.86 eV above zero), but as these calculated lifetimes are more than an order of magnitude too large it is safe to accept that the natural width of a transitioning state near threshold will still far exceed the spacing between the states.Hence, the motions in the thermal reaction range can be assumed to be chaotic, and a classical simulation should mimic the true behaviour.
Were this not to be so, Baer and Hase 17 consider a limiting case where the phase space is divided into two regions, "one part consisting of chaotic trajectories which can decompose and the other of quasi-periodic trajectories which are trapped in the reactant phase".They then introduce the ratio f ch ¼ r ch (E)/ r(E), being the fraction of phase space that supports chaotic motion, and denes the amount by which the computed classical rate would fall below the true value.

Computational procedure
For a CH 3 NC or CH 3 CN molecule in its normal internuclear conguration, each of the 6 atoms is assigned a set of random momenta in the x, y, and z directions; any unwanted translation and rotation are removed, and the motions scaled to the target energy value.Then, for each initial energy distribution, aer a simple fourth-order Runge-Kutta start-up procedure, the trajectory is followed by a h-order Adams-Mouton predictorcorrector.For each trajectory, a standard time step length of 0.1 fs was used, and arithmetic sums of all N values of internuclear displacements i.e.

RðtÞ¼
X N i¼1 È rðtÞÀ r e É i were tabulated at dt ¼ 2.5 fs intervals; † for some applications, a vector P(t) of the instantaneous potential energy values was also saved.Throughout this study, the minimum ensemble size was 1000 trajectories.These methods have been described extensively in previous publications, 12,13 and all algorithms used are described in detail elsewhere. 18n the preceding paper, 13 I speculated that a possible surrogate for f ch can be obtained from the spectral properties of such classical trajectories.The Fourier transform of the vector R(t) yields an overall spectrum of the collective vibrations.The spectral resolution is determined by the vector length, and with 32 768 (2 15 ) terms yields a resolution of about 6.5 cm À1 , thus requiring trajectory lengths of 81.92 ps.The area under each spectral curve is a measure of the total spectral intensity and rises linearly with increasing energy above the classical minimum. 13The postulate then is that the fraction of the total intensity lying below the lowest vibrational frequency, 263 cm À1 for CH 3 NC, 19 gives a "measure" of the degree of chaos in the system (and for simplicity is taken to be the area between 0 and 250 cm À1 , expressed as a percentage, and called f c ).That is, if the area under a regular spectrum is proportional to the total vibrational energy, then any energy that is not "seen" by the Fourier transform must be chaotic.As the degree of chaos increases, and the motions degenerate into white noise, for which the Fourier transform is at, the spectral features shrink and the missing energy migrates to a singularity at zero frequency; the properties of this "spike" are that the greater the number of points in the given frequency range, the narrower and taller this feature becomes, tending to innity in the limit. 20his behaviour is illustrated in Fig. 1: part (a) shows a typical spectrum for a trajectory that does not isomerise and by comparison with Fig. 3 of the earlier paper, 13 which shows the transforms of the 3 C-H or the 3 H-H distances, these features are related to the motions within the CH 3 group and are still partially regular; part (b) shows another trajectory, at the same energy, where the spectrum is at.Instead of terminating this trajectory at the rst crossing, it was allowed to run on to 81.92 ps, whence it made 37 crossings, sharing its time between CH 3 NC and CH 3 CN; the crossings occur at random potential energy values between 2.9 and 5.4 eV, so that each newly formed isomer evolves from a totally different set of positions and momenta, giving quasi-periodic snippets strung together in random fashion, ensuring the illusion of complete chaos.Also, as dots, is the same spectrum expanded by a factor of 5000, where the remnants of the CH 3 NC and CH 3 CN spectra still persist.This is a much stronger demonstration of a chaotic spectrum than found previously. 13ig. 2 shows the behaviour of these f c values for CH 3 NC from an energy of 2 eV above the global minimum (in the classical region, just below the zpe) up to 7.3 eV (zpe + 5.1 eV), at which point, only 20 out of 1500 trajectories survive for 81.92 ps; this is an exceedingly high energy, being 115 kcal mol À1 above the zpe whereas the activation energy is only 38.2 kcal mol À1 . 16As these trajectories are only partially chaotic, there is no reason for these chaos fractions, f c , at any energy to coincide since each trajectory starts from a unique energy distribution and retains some memory of it.Fig. 2 also shows the corresponding curve for CH 3 CN, starting at its zero-point energy; up until a total energy of 5 eV, the pattern is very similar to that for CH 3 NC, but with a downward shi of about 1 eV, the difference in energy of the two systems.Beyond that, there a sharp turn upwards, but it was not possible to collect a set of comparable values at the highest energy as out of 5000 trajectories only one lived for 81.92 ps.

The standard method
The trajectory observations are very similar to those reported previously, 21 but then trajectory lengths were restricted (by a CPU speed of 1 MHz) to 3 ps, whereas now, 1 ns is readily accessible.Much of the analysis in this Section uses trajectories for energies of 4.0 and 5.0 eV above the CH 3 NC potential minimum because, as just noted, too few trajectories survive for the required time at higher energies; conversely, as the energy is reduced, fewer trajectories react, and eventually, none at all.
As we, 21 and others, 22 found, there is a sharp drop within the rst ps in the number of molecules remaining, due to some possessing excess energy in the critical directions reacting before equilibration could occur.Then (at a similar energy), a long period, up to over 200 ps (instead of the original 3 ps) 21 when the logarithm of the number remaining is closely linear with respect to time; aer that, the points become rather scattered, as those few of the original 1000 now le no longer comprise an acceptable statistical sample.
Although it may seem plausible, the lifetime of a trajectory at any given energy does not correlate with its degree of chaos, f c : for example, at 4.0 eV above the CH 3 NC potential minimum, the average f c for trajectories with lifetimes between 89 to 222 ps is 6.6% and for 275 to 485 ps, it is 6.5%.
It should also be noted that the spectrum shown in Fig. 1(a) and all those used in the calculation of f c are different from the composite spectra shown by Sewell and Thompson: 23 those are the sum of the Fourier transforms of each bond distance, whereas these are the Fourier transform of the sums of all interatom distances; these are not the same although, oen, they can look similar if long-range interactions are unimportant.
Thus, we have a rst-order decay process, giving the appearance of a conventional unimolecular reaction, 17,24 although none of these single-pass trajectories (as in Fig. 1(a) and similar ones at all lower energies) is anywhere close to being chaotic; otherwise, these spectra would be at.As seen in Fig. 2, the prospective measure, f c , of the degree of chaos rises from near zero at low energy to between 5% near the activation region and 10-15% at unrealistically high degrees of excitation.
First, we might consider this as an extreme limiting case of "intrinsic non-RRKM behaviour", as depicted in Fig. 1(e) of Bunker and Hase, 24 and discussed in more detail by Baer and Hase later. 17Aer the initial loss of specically excited molecules during the rst ps, we are le with two non-interacting regions of phase space, one whose properties are unknown and does not decay, and one which is chaotic and decays in the expected manner; this could account for the very low rate constant.
In fact, we nd that the whole CH 3 NC ensemble reacts before 385 ps at 5.0 eV, or before 485 ps at 4.0 eV above its minimum.Now add to Fig. 1(e) a slow leakage between the chaotic and the quasi-periodic regions, which Bunker and Hase actually show in their Fig.4(c), 24 but never discuss.From the manner in which the initial energy is allocated to the molecules, very crudely,  about two-thirds of the total momentum will reside in the CH 3 moiety.During the rst ps, 21,22 every molecule assigned sufficient N]C angular motion will react, leaving a sub-ensemble of molecules with mainly C-N]C vibration and a modicum of C-N]C bending motions which are isolated from the motions inside the CH 3 group.The excess energy in the methyl group then leaks rather slowly (on a roughly $100 ps time scale) by a rst order process whence the N]C group attains sufficient angular motion to invert almost instantly.Because of the large difference between the rates of these two processes, the population versus time plot in Fig. 1(e) collapses into bimodal form, a sharp initial drop followed by a rst order decay.While originally these calculations were assumed to yield a reaction rate, 21 all that was being observed was a spurious internal relaxation!
The disparity between the calculated and the measured rates is shown in Fig. 3, the solid line being an estimate of the specic rate function, k(E), derived as the Slater-Forst inversion of the experimental Arrhenius plot, 25 and the bottom set of discrete points being the calculated results; the open circles are the new and the small solid squares the earlier 10 results with, respectively, 1000-and 100-trajectory ensembles.

Further analysis
It is not clear what causes this difficulty but the question should be resolved for this particular reaction once an ab initio potential energy surface becomes available.But, for now, it is easy to rule out the computational method, or any lack of numerical precision, as the cause of this huge degree of disagreement.
Throughout the history of such calculations, it has been known that the precision of the computed trajectories has seemed almost irrelevant.Hase and Buckowski 26 used the reversal of trajectories to test whether they would return to the original starting conguration, which they did not, and we have demonstrated many examples of vastly different individual trajectory lifetimes with the same algorithm on different machines, or different operating systems on the same machine; 12,27 despite such variation in lifetime for the same trajectory, sometimes in excess of a factor of 200, the ensemble lifetime averages agreed well, even for quite small ensembles. 26,27What is important in these calculations is not precision in following the trajectory itself, only that the trajectory should be able to go anywhere within the allowable phase space as long as the energy and angular momentum are conserved; this is where precision is vital. 12o the aw must lie with the potential energy surface, despite its apparent well-behaved properties.For example, it may be that although this surface conforms to all the parameters, both at the normal structure and near the asymptotic limits, there may be some intermediate impediment that tends to impose "intrinsic non-RRKM behaviour" on the computed solution.It may be that the surface enforces too much quasi-periodic motion within the CH 3 group, and thus an insufficient degree of chaos; this would be consistent with f c values of the order 3-5% in Fig. 2 and the decit of a factor of around 20 in the calculated k(E) values over the actual experimental energy range, in Fig. 3.
One proposal for increasing the degree of chaos in a trajectory was to make deliberate discontinuities in the integration paths, chosen so as not to compromise the conservation of energy.A simple way is to reverse the momenta of a few chosen particles at regular time intervals, e.g.0.1 ps, 10 but this only conserves the energy, and not the angular momentum unless it is zero; however, methods for conserving both have been explored previously. 28In this way, the computed k(E) value for the CH 3 NC / CH 3 CN reaction was raised by a factor of about 15 near reaction threshold, and of about 5 at very high energy. 10n that study, only reversals of peripheral atoms were used, the N atom and either one, or all three, of the H atoms, and the latter was by far the more effective in accelerating the reaction rate.Here, not every permutation of possible groupings was examined, but it appears that the original 4-atom switch (i.e.N and all H atoms) is the most effective.As an aside, reversal of all 6 atoms results in no reaction, as the system meanders about the initial reversal point, and only dris away slowly due to rounding errors in the 16th decimal place; also, reversals of 5 atoms each time, or of two groups of 3 atoms alternately, are less effective.
The results of this reversal process are also shown in Fig. 3 as open squares: at high energy, they agree well with the older ones, 10 shown as solid triangles, but with a steadily increasing difference at lower energy.‡ In contrast, the mean lifetimes for the reverse reaction CH 3 CN / CH 3 NC at the two highest energies are unaffected by  10 ‡ This probably arises because the new reversal process is programmed with fewer instructions, thereby minimising the growth of computational chaos; the new results are preferable as the older ones place the low-energy rate well above the experimental estimate.the reversal procedure; at lower energies there is a small acceleration, suggesting that impediments to reaction occur mainly on the CH 3 NC side of the potential energy surface.This observation also helps to resolve the microscopic reversibility problem: the approximate state densities at 41 150 cm À1 above the global minimum are 12 r CN ¼ 2.4 Â 10 7 and r NC ¼ 6.7 Â 10 6 states per cm À1 respectively, then r CN /r NC x 3.6 which should be equal to s CN /s * NC ¼ 8.1 (ps)/3.3(ps) x 2.5; here, the asterisk denotes the reversal procedure and given the uncertainties of at least AE5% in the 1000-ensemble lifetimes and of those in the state densities, this agreement is promising.
Finally, if the original data (open circles in Fig. 3) are scaled by the corresponding value of f c as suggested by Baer and Hase, 17 we get s Â NC ¼ s NC /f c , depicted as crosses in Fig. 3.The s Â NC values track the s * NC points fairly closely, giving s CN /s Â NC ¼ 8.1 (ps)/2.4(ps) x 3.4, but detailed analysis is not yet warranted because of the crudity of the state-density estimates; at this point, an exact count would be needed, together with a detailed analysis of the anharmonicities of all the vibrations in both molecules.

Summary of reversal patterns
In this Section, we need to calculate several thousand trajectories to illustrate a single item: hence, most of the following discussion is conned to higher energies, principally 7.0 eV above the CH 3 CN minimum (6.0 eV above the CH 3 NC minimum).
(i) Without reversals: Fig. 4 shows the average distribution of potential energy for both CH 3 NC and CH 3 CN at a total energy E ¼ 7.0 eV.The upper histograms give plots of the vector P(t), with the instantaneous potential energy values sorted into 0.1 eV bins.As we have shown previously, 11 the potential energy distribution for a collection of non-interacting harmonic oscillators can be represented by and it was found to be approximately true for the CH 3 NC molecule, with a value of m x 4.3, even though some of the oscillators were either degenerate and/or anharmonic.In Fig. 4(a), the histogram is overlaid by this sine function with m ¼ 4.3, which is quite a good t; the same is true for CH 3 CN in Fig. 4(b), but the shapes appear different because in (a), P is scaled to lie between 1.0 eV and 7.0 eV, and in (b), between zero and 7.0 eV.In fact, this sin m (pP) shape with m ¼ 4.3 ts all potential distribution histograms within 4.3 # E # 7.3 eV for both molecules.
Also, in both of these diagrams, downward-pointing zpe markers are placed at 1.2 eV above their respective minima, showing that either molecule spends very little time in the forbidden region below the zero-point energy; however, this fraction increases rapidly as the total energy is reduced, but once the value of m for the potential energy surface being studied is known, the fraction is easily estimated with the above formula.
The lower histogram in Fig. 4(a) is that of the potential energies at which 2000 trajectories crossed from CH 3 NC to CH 3 CN; the vertical markers indicate the range over which P is scaled, the height of the transition barrier (2.8 eV), and the total energy E (7.0 eV).In Fig. 4(b), the corresponding histogram is that for the opposite reaction, and the bottom curve is a copy (slightly displaced) of that from Fig. 4(a) showing that these two proles are the same.This observation implies, importantly, that if neither set of molecules has any impediment in accessing the transition region, microscopic reversibility will be obeyed.
The most puzzling aspect of these two reaction potential energy proles is that they, too, conform to a sin m (pP) shape, but the best t is for m x 12, which happens to be the number of oscillators in each system.However, for n uncoupled oscillators, 11 m x (0.8n À 0.7), n $ 4 so this coincidence could be an artefact of the potential surface in the vicinity of the transition point.It should be noted that for the reaction NCNC / NCCN, using an ab initio potential surface, the overall and the reactive potential energy proles gave approximately the same value of m, 11 and not widely different values as we see here.
(ii) With reversals: Fig. 5(a) and (b) show corresponding graphs for trajectory ensembles in which the motions of the N atom and all three H atoms were reversed every 0.1 ps.Now, the potential energy distribution functions exhibit the same values of m (4.3 and 12) as do the unperturbed histograms in Fig. 4, but all the curves are shrunk downwards in energy by $0.5 eV compared with the sin m (pP) functions copied from Fig. 4. It appears that, somehow, this procedure prevents access to the topmost energy range.
More detailed examination reveals, in fact, that the reversal procedure does not alter the degree of chaos in these systems.If the degree of vibrational mixing, and therefore chaos, were to be increased by the reversal process, then the value of m tting the overall potential energy distribution would be smaller in Fig. 5 than in Fig. 4. 11 Moreover, the values of f c found in the Fourier transform spectra of R(t) remain within the same bounds as in Fig. 2, with or without reversals.Clearly, the introduction of additional computational chaos does not increase the degree of vibrational chaos; it seems that it must simply increase the chance of the trajectory nding itself in the vicinity of the transition region.However, how it reduces the mean potential energy of the system, and of the reacting molecules, is not obvious.
Fig. 5(a) includes an extra variation, that of random reversal spacings in the range 0.1-0.2ps, to eliminate the possibility of a spectral feature related to the reversal frequency, but none was found.Both the overall potential energy and the reaction energy prole (bottom line) remained the same, but the mean lifetime was about 15% longer, possibly because the mean time between reversals was now 0.15 ps instead of 0.1 ps, with a concomitant reduction in computational chaos.
A number of mainly negative conclusions can be deduced from these two pairs of diagrams.For both CH 3 NC and CH 3 CN, the reversal procedure causes the overall PE and the reactive PE curves to squeeze down to lower energy by about 0.5 eV; in case Fig. 5(a), the mean lifetime at total energy 7.0 eV is reduced from 20.8 to 3.7 ps, but in case Fig. 5(b) only from 11.6 to 11.5 ps.As the relative shapes and positions of the normal and reactive PE proles hardly change, the large lifetime differences do not depend directly on these factors.

Conclusions
Accepting that, for the foreseeable future, only semi-empirical potential energy surfaces will be available for estimating reaction rate constants of many interesting molecules, this study begins to explore methods for validating such surfaces.The example chosen is that of a very simple reaction, for which it has long been known that standard classical trajectory methods yield grossly incorrect rate constants.Only two facets of the problem are addressed, whether energy can ow freely between different parts of the molecule, and the energy prole of those molecules reacting relative to the whole energy distribution prole.
Despite the simplicity of the two avenues explored here, some aspects of our understanding have been improved, and a number of new ones have been observed, as follows: (a) It is shown that the proposed heuristic measure, f c , of the degree of chaos behaves appropriately as the total chaos limit is approached -Fig.1(b).
(b) With an improvement in the earlier 13 calculation of f c , the scatter of values for CH 3 NC is reduced, and the resulting degree of chaos rises almost linearly with increasing energy; however, the corresponding curve for CH 3 CN tends to rise much faster at higher energy and may be approaching the 100% limit -Fig.2.
(c) For a given energy, there is no correlation between the lifetime and the degree of chaos, f c .
(d) The fact that all trajectories eventually react means that this particular PE surface does not impose strict "intrinsic non-RRKM behaviour" but a modied version of it.
(e) Assuming the approximate equivalence of f c with the f ch of Baer and Hase, 17 scaling the computed lifetimes with the inverse of f c yields rate coefficients fairly close to the experimentally derived values -Fig.3.
(f) It also appears to resolve the problem with microscopic reversibility.
(g) A deliberate attempt to amplify the degree of computational chaos by periodic reversal of momenta achieves approximately the same result -Fig.3.
(h) Application of the same reversal procedure to CH 3 CN trajectories causes little or no acceleration at high energies, consistent with the high-energy behaviour of f c in Fig. 2, meaning that the reverse reaction rate is close to the required value; this implies that reaction impediment is very much weaker on the CH 3 CN side of the potential barrier than it is for the forward reaction.
(i) For all total energy values studied, and for both CH 3 NC and CH 3 CN, the distribution of potential energy followed a sin m shape, with m $ 4.3; likewise, the potential energy of reacting molecules also followed a sin m shape, but with m $ 12 -Fig.4.
(j) The fact that m $ 12 for both forward and reverse reactions proves that microscopic reversibility must hold if access to the transitional congurations is equal on both sides of the barrier.
(k) Inclusion of momentum reversals every 0.1 ps at 7.0 eV above the global minimum moves the sin m curves downward by about 0.5 eV but does not change the values of m; hence, computational chaos does not affect dynamical chaos -Fig.5 and discussion thereof.
(l) Downward movement of the mean potential energy towards the saddle point accelerates the reaction: hence the defect(s) in the PE surface must be in the mid-energy range, but mainly on the forward side.
(m) Choice of random intervals between the momentum reversals shows that no spurious spectral features arise, but that the longer the average interval, the longer the average lifetime.

Discussion
It is shown decisively that for the test reaction, although the potential energy surface was carefully craed, there is a strong impediment to the ow of energy for CH 3 NC molecules to access the transition state, but a much weaker one for the CH 3 CN molecule to do so.Once these impediments are removed, the equivalence of the reactant energy proles for forward and reverse processes guarantees that microscopic reversibility will ensue.This will be possible once an ab initio potential energy surface becomes available, by showing which volumes of conguration space are being excluded; recent machine-learning computer methods 29 could be invaluable in the remediation of defects in this PE surface.
The adjustments needed in this case could also give clues for the improvement of other similar empirical surfaces where ab initio data are not available, for example the dissociation of the peroxy radical CH 3 -O-O-N-O / CH 3 O + NO 2 studied by Stimac and Barker: 30 they also found the rate a factor of 20 too small but, as well, the ensemble decay resembled that of a triatomic molecule. 31Of course, this factor of 20 could just be a coincidence as the potential energy surface is far more tangled, and there might be signicant complications from "roaming" effects, such are found in the dissociation of other nitrocompounds. 32We should note that for a dissociation reaction, the comparison of forward and reverse reactions is more difficult than for a simple isomerisation process. 33s already suggested, the scaling of computed k(E) values with 1/f c appears attractive, but it is only by comparison with results when using a true ab initio potential surface that it can be shown to be valid.On the other hand, the use of reversals to impose massive amounts of computational chaos, articially planting trajectories on inaccessible reaction paths seems to provide a crude, but extremely powerful, qualitative mechanism for detecting bottlenecks in the activation process.How it distorts the potential energy proles in Fig. 5, and thereby reduces the trajectory lifetimes is very puzzling; however, this reversal procedure does not represent any real physical process, so the results may remain as an intellectual curiosity for now.
Possibly more instructive are the two kinds of potential energy plots.For the overall PE distribution, it is a simple matter to save the values of the energy at a large set of times for any non-reacting trajectory and sort the values into bins.However, for the reactive PE proles, a very minimum of 2000 trajectories is needed to generate a useful histogram.These proles have three important uses: (i) the forward and reverse proles must agree; (ii) the existence of only one simple sin m feature may help to guard against an unexpected alternative reaction channel; and (iii) that the values of m for the overall and the reactive proles should not differ wildly, as the do here.Item (iii) is a matter of speculation, based on the observation of the similarity of the two values of m found in the only known calculation of this kind using an ab initio potential energy surface. 11xtension to larger molecules adds some complications.In addition to the increased complexity of the potential function, the number of inter-atomic distances needed to calculate R(t) increases rapidly, e.g. for a 10-atom system there are 45 distances, but present results suggest that only nearest and next-nearest neighbour interactions make any signicant contribution to the spectrum itself, or to f c ; however there could be large but compact molecules, like neopentane or ditert-butyl peroxide, where this may not be the case.
Finally, the zero-point energy problem has to be addressed.38,39 Recently, the rst approach has been exploited beautifully to describe the internal motions of small water clusters. 40,41However, this degree of sophistication for large organic molecules, projecting the motions into normal modes could be prohibitive, but a simple trial-and-error momentum reversal of one or a few atomic motions to prevent crossing below the zpe could probably be found.
The second approach has been shown to be exact in the RRKM limit, 39 and the time fraction can be followed easily during the trajectory calculation, estimated by state-counting methods, 10,11 or simply from the potential energy distributions as shown in Fig. 4.
In either case, if it is thought that tunnelling may be an issue, trajectories may penetrate the potential barrier below its maximum, and protocols for decay of the associated wave packet inside the classically forbidden barrier region need to be established. 10o summarise: momentum reversals are a useful tool for detecting some energy redistribution aws in potential energy surfaces, but it is unclear how this works; attempts to quantify the degree of chaos show some interesting aspects but, for now, the results are tentative; potential energy distributions, as shown in Fig. 4, are especially informative for isomerisation reactions.

Fig. 2
Fig. 2 Percentage of spectral intensity residing below 250 cm À1 versus energy above potential minimum; the vertical bars on each point represent the standard deviations from a random selection of 20 trajectories.The position of the barrier for both reactions, and for the zpe of CH 3 NC are shown by vertical lines; the graph for CH 3 CN starts at the zpe.

Fig. 3
Fig. 3 Computed specific rate function for CH 3 NC / CH 3 CN versus energy above zpe.Line -Slater-Forst inversion of experimental Arrhenius rate law; circlesplain calculation; open squares -4-atom reversals at 0.1 ps intervals; crossesopen circles scaled by f c .Solid squares and triangles are, respectively, plain and 4-atom reversal results from earlier calculations.10

Fig. 4
Fig. 4 Potential energy distributions for a total energy of 7.0 eV.The large top plots are histograms of the 32 768 energy values P(t); the small lower plots are histograms of the potential energy crossing values for 2000 trajectories each, and the lowest one in (b) is copied from panel (a) to demonstrate their equivalence.The smooth curves are fitted sin m (pP) functions with m ¼ 4.3 or m ¼ 12, as described in the text.

Fig. 5
Fig. 5 Effect of equally-spaced reversals on the potential energy distributions, P(t), for a total energy of 7.0 eV.The smooth curves are copied from Fig. 4(a) and (b) respectively, and the two pairs of histograms show a shift downwards of $0.5 eV in each case.The lowest histogram in panel (a) shows the crossing energies of 2000 trajectories with randomly-spaced reversals; as in Fig. 4, all fitted sin m (pP) curves have either m ¼ 4.3 or m ¼ 12.

Table 1
Mean lifetimes (s NC ) of classical trajectories for isomerisation of a (J ¼ 0) methyl isocyanide ensemble at 3 different excitation energies