Arup Chakraborty Research Group

A computational framework capable of predicting reaction pathways and rate constants would be extremely useful in many applications.  Reactions typically involve bond making and bond breaking which is not easily described by model potentials.  This difficulty necessitates a method based on ab initio potential energy calculations.  Ab initio molecular dynamics (MD) simulations (such as Car-Parinello) are impractical for studying reactions because reactions are rare events, i.e. one reaction event is separated from the next by a time much longer than the duration of a feasible MD simulation.  It is equally impractical to exhaustively map out an entire potential energy surface even for systems with as few as four atoms.¹  What is required then is an efficient way of locating and describing critical portions of the potential energy surface, obtaining reaction coordinates, and finally computing rate constants. 

Ideally, a calculation of this sort should start from a reactant state² and calculate – without relying on intuition – the relevant transition states, reaction pathways, transition state theory (TST) rate constants, and non-TST corrections.  There are six overall steps in this procedure.  Each step is described below with relevant references and a pictorial example where possible.

1)         Locate saddle points by an automated Cerjan-Miller type search procedure³

 2)    Determine the reaction pathway associated with a saddle point^4,5,6.   This involves computing the steepest descent reaction path and the transverse eigenmodes and frequencies along the steepest descent path.  From this information, the reaction coordinate and an approximate potential energy surface are identified.  The figure below at left shows the true Muller surface (with stationary points in red).  The middle figure is the resulting approximation to the Muller surface, and the right figure shows difference between the true and approximate surfaces where the difference is within 1% of the saddle point to reactant energy barrier. 

3)      Perform additional saddle point searches, excluding the pathways already found so as to locate multiple reaction pathways available

to the reactants.  The figure at right illustrates the

construction of an umbrella potential to block the

pathway we located above on the Muller surface. 

The umbrella potential is shown above the reaction

pathway.

4)      Calculate a variational TST rate constant for each pathway⁷.  In the figures below, free energy barriers for cyclohexane boat – chair interconversion have been computed using the method described.  The potential energy surface that gave the results at left is an empirical potential energy surface^8,9,10 and the results at right are from Hartree-Fock calculations¹¹. The variational TST rate constant is (kT/h)exp(-DF_max/kT).  In both cases the red dots are the potential energy along the steepest descent path.

5)      Calculate a reaction path Hamiltonian (RPH) for each pathway¹.  In the RPH, the reaction coordinate is naturally coupled to the transverse (bath) modes due to curvature of the reaction pathway and twisting of the transverse modes about the pathway.  The coupling terms in the Hamiltonian require extremely accurate description of the steepest descent path. 

6)      Calculate corrections for non-TST effects such as re-crossing and tunneling.^12,13,14,15  This requires dynamic simulations of the barrier crossing events which can be done using Hamilton’s equations of motion. 

 This method has been applied to several small systems including H-transfer in malonaldehyde, ammonia umbrella inversion, and HCN rearrangement to HNC.  Currently this method is being applied to study proton mobility in H-ZSM5.  I am coadvised by Professor Alex T. Bell (www.cchem.berkeley.edu/~atbgrp/) and Professor Arup K. Chakraborty.  All electronic structure calculations in this work are done using QChem 2.0 (www.q-chem.com).  This project is shared by Dr. Wan Zhen Liang in the Chakraborty group and myself.

 References

1                     W.H. Miller, N.C. Handy, J.E. Adams, J. Chem. Phys., 72, 99 (1980).

2                     www.q-chem.com an example of a commercial code with a built in optimization algorithm for finding reactant configurations.

3                     C.J. Cerjan, W.H. Miller, J. Chem. Phys., 75, 1280 (1981).

4                     E.M.Sevick, A.T. Bell, D.N. Theodorou, J. Chem. Phys. 98 3196 (1993).

5                     M. Page, C. Doubleday, J.W. McIver, J. Chem. Phys. 93 5634 (1990).

6                     M. Page, J.W. McIver, J. Chem. Phys. 88 922 (1988).

7                     D.G. Truhlar, W. L. Hase, J.T Hynes, J. Phys. Chem. 87, 2664 (1983).

8                     H.M.Pickett, H.L.Strauss., J. Amer. Chem. Soc. 92 7281 (1970).

9                     R.A. Kuharski, D.Chandler, J.A. Montgomery, F. Rabii, S.J. Singer, 92 3261 (1988).

10                 M.A.Wilson, D Chandler, Chem. Phys. 149 11 (1990).

11                 www.q-chem.com

12                 D. Chandler, U.C.Berkeley Stat. Mech. Course Notes, Odin, Berkeley (2000).

13                 D.Chandler, Introduction to Modern Statistical Mechanics, Oxford University Press, New York (1987).

14                 W.H.Miller, J. Phys. Chem. A 102 793 (1998).

15                 P. Hanggi, P. Talkner, M. Borkovec, Rev. Mod. Phys. 62 251 (1990).