Rice–Rampsberger–Kassel–Marcus (RRKM) theory calculates an energy-dependent microcanonical unimolecular rate constant for a chemical reaction from a sum and density of vibrational quantum states. This article demonstrates how to program the Beyer–Swinehart direct count of the sum and density of states for harmonic oscillators, as well as the Stein–Rabinovitch extension for anharmonic oscillators. Microcanonical rate constants are calculated for the decomposition of vinyl cyanide () into , and as an example.

1. Introduction

The essential framework for our understanding of unimolecular reaction rates was developed largely in the first half of the twentieth century [1]. The Lindemann–Hinshelwood mechanism, proposed in 1922 by Lindemann and later expanded upon by Hinshelwood [2], posited that the gas phase unimolecular reactions were initiated by bimolecular collisions. This allowed the rate of the reaction to be calculated as a function of internal energy, regardless of the method of activation. Later that decade, Rice and Ramsperger and, independently, Kassel developed an improvement on the model known as the Rice–Ramsperger–Kassel (RRK) model [3, 4]. The RRK model viewed the molecule as a system of identical harmonic oscillators and introduced the idea of activation energy: that sufficient energy must be deposited into specific modes of motion in order for the reaction to occur. In 1952, Rice and Marcus extended this into what is now known as Rice–Rampsberger–Kassel–Marcus theory (RRKM), which incorporates a more complete quantum mechanical description of the molecule and utilizes the concept of the transition state, a specific conformation that the molecule must adopt in order for the reaction to proceed [5]. The RRKM model remains the prevailing scientific description of unimolecular chemical kinetics to this day.

Chemical transformations require bonds to break in the reactants and new bonds to form as products are generated. The potential energy associated with these transformations defines an -dimensional potential energy surface (PES), where is the number of atoms [6]. The reaction pathway is a one-dimensional cross section of the PES composed of the precursor, intermediates, transition states and products of a particular reaction [6]. Reactants and intermediates of a reaction correspond to minima, whereas transition states are located at maxima along the reaction pathway.

The rate of the reaction depends on the probability of the molecule adopting the conformations along the reaction pathway, which, in turn, stems from the probability that a sufficient amount of energy is partitioned into the necessary modes of motion. The RRKM formalism assumes that: (1) once the molecule adopts the transition state orientation, the reaction proceeds to products; and (2) the internal vibrational energy redistribution (IVR) is fast with respect to the timescale of the reaction. Under these constraints, the RRKM microcanonical rate constant, , is given by

(1)

where is the internal energy of the system, is the activation energy for the reaction, is the sum of states of the transition state from to , is the density of states of the precursor at energy , is Planck’s constant, and is the degeneracy of the reaction pathway [3, 4, 5]. A state is defined as any unique vibrational configuration that determines the internal energy of the molecule. The sum of states, , is the total number of states within a specific energy range [6]. The density of states is the number of states per energy level. Equivalently, the density of states is the derivative of the sum of states with respect to energy [6]. The density is expressed in units of inverse energy
(), while the sum of states is a dimensionless quantity. Thus, the microcanonical rate constant of equation (1) has units .

There are vibrational frequencies for nonlinear polyatomic molecules and vibrational frequencies for the corresponding transition state. It is the values of these frequencies, which typically range from , that dictate how the internal energy is distributed in the molecule. The calculation for the sum and density of states depends on these frequencies, which can be measured or looked up in standard data tables. However, for transition states, these values must be estimated or, more likely, calculated using theoretical packages.

This article is divided into the following sections: In Section 2, the sum and density of states for vinyl cyanide () are determined through a direct count method analogous to that developed by Beyer and Swinehart [7]. In Section 3, the method is extended to anharmonic oscillators using the Stein–Rabinovitch algorithm [8]. In Section 4, the RRKM microcanonical rate constants are calculated for the unimolecular dissociation of vinyl cyanide and are used to predict various dynamic properties of the chemical reaction.

2. Direct Count for Harmonic Oscillators: The Beyer–Swinehart Algorithm

Beyer and Swinehart (BS) developed a surprisingly simple yet computationally intensive direct count method to calculate the sum and density of states at defined energies of a harmonic oscillator [7]. The advantage of conducting such calculations in Mathematica is that the user can leverage high-quality numerical ODE solvers and interactive graphical features to visualize dynamic properties of the chemical reaction. This is demonstrated in Section 4, where the temporal development of nine different species is plotted and automatically updated as the user sets the internal energy of the system with a slider. All of this can be done in real time after some initial computations that are exact with Mathematica.

There are a number of equivalent ways to implement the BS direct count of states. To ease the exposition, we first illustrate the task with a loop. The sum and density of states are determined by the vibrational frequencies of the species of interest, which we define in the list . In more complex situations, such as biomolecules or reactions at very high energies, to lessen the computational burden it may be necessary to divide the value of the total energy into larger packets and modify the vibrational frequencies to be multiples of the new packet size [8]. In this case, the packet size is 1 and the vibrational frequencies for vinyl cyanide are taken from reference [9].

The final result for the sum or density of states is a list with elements. For the sum of states, this list is initialized with all values set to 1. The sum of states can then be calculated using a nested (or doubly indexed) loop.

Because the density of states is the derivative of the sum of states, the density of states can also be determined from the sum of states table by setting the element of the density of states to the difference between the and elements of the sum of states table. Alternatively and equivalently, this can be done to the initial sum of states table to derive the initial density of states table, which is a list with elements with the initial element set to 1 and all other elements set to 0. This can be conveniently done with .

A visual comparison of the sum and density of states, calculated for vinyl cyanide (), shows the behavior of the two quantities over the energy range. As the general trend of the graphic is clear, we subsample the lists of values to avoid plotting millions of points, needlessly increasing file size.

Thus, when vinyl cyanide has an internal energy of , the sum and density of quantum states are approximately and .

Notice that the loops for computing the sum and density of states perform the same basic process—the BS direct count method—on different initial lists. As we proceed, adopting a functional programming approach eases the exposition and vastly improves code clarity. We do this by introducing two new functions, and . The first updates the lists of sums and densities for a single vibrational frequency, and the second uses to iteratively update for each of the vibrational frequencies in the table.

We then use these to create the two functions and , which take a frequency list and total energy as arguments and do just as their names imply. The last argument to these functions, , is used in Section 4.

For completeness, we recompute the sum and density of states using these functions to illustrate how that can be done. It is readily verified that these compute the same quantities.

3. Direct Count for Anharmonic Oscillators: The Stein–Rabinovitch Extension

The harmonic oscillator model assumes that the difference between successive vibrational energy levels remains constant. In real oscillators, this difference decreases as the vibrational excitation increases. These anharmonic effects increase both the sum and density of states and thus often cancel when taking the ratio to determine the RRKM rate constant in equation (1). However, there are some cases where this effect does not completely cancel or where highly resolved measurements necessitate a more accurate description of the vibrational energy.

The Stein–Rabinovitch (SR) algorithm extends the BS method to incorporate the effects of anharmonicity into the RRKM rate constant [8]. Following [10], the anharmonic effects are incorporated into the energy level expression through the equation

(2)

where is the energy of the vibrational level, is the vibrational frequency, and is the anharmonicity constant. The subscript stands for equilibrium. It is chosen because this would be the vibrational frequency if the oscillator vibrated harmonically about its equilibrium position (the minimum in the well). The SR extension counts each energy level of the anharmonic model to calculate the sum and density of states in a way that is functionally similar to the BS method; however, the input table is modified to include the anharmonicity constant for each vibrational mode. For this example, we have organized the anharmonicity constants, , into the table . Using the vibrational frequencies of vinyl cyanide from [9], and assigning 0.01 as the anharmonicity constant for each vibration, the energy for each level of each mode can be calculated.

With these we can calculate the sum and density of states, initializing the lists in the same manner as the BS method (a list of 1s for the sum of states, and a 1 followed by 0s for the density of states). The process is functionally similar to that used to calculate the sum and density of states for a set of harmonic oscillators. In the BS method, the elements of a single table are modified with each iteration of the inner loop. By contrast, in the SR extension, the elements of the table are updated through separable computations combined after each complete loop for each oscillator.

The implementations can be verified by replicating the calculations in the appendix of [8], which shows the calculation of the density of states up to 10 units of energy for a model system with two oscillators. The energies of the first oscillator are and the energies of the second are .

With these, we can define SR analogs to and .

We can then use them to compute the SR sum and density of states.

Here are two plots comparing the harmonic and anharmonic models for the density and sum of states calculated for vinyl cyanide.

4. Dissociation of Vinyl Cyanide

As a more substantive example, we calculate rate constants for the decomposition of vinyl cyanide () into hydrogen cyanide (), hydrogen isocyanide () and acetylene (), compare them to published values, and interactively visualize the temporal dependence of the various species. In [9], Homayoon and colleagues considered the following kinetic model that involved seven different paths, three different intermediates, 13 rate constants and three different products:

Each step of each path passes through a transition state to yield either an intermediate or a final product. The following differential rate equations characterize the temporal dependence of the amounts of each species:

The determination of the sum and density of states requires the vibrational frequencies for each species located at minima and maxima along the reaction pathway. For the decomposition reaction of vinyl cyanide, there are six minima (the geometry of along with each intermediate) and 11 maxima (each unique transition state). These 17 sets of frequencies are taken from reference [9] and are provided in Table 1. In the following definitions, we have chosen the same notation as in reference [9]. TS indicates a transition state, the particular path is indicated by a Roman numeral, and the minimum energy of each species relative to the minimum energy of vinyl cyanide has the prefix re.

We now turn to the calculation of the sum and density of states and then use the values to calculate each rate constant. We begin by defining the vibrational constants.

Because the internal energy is defined relative to the well depth of vinyl cyanide, the sums and densities of states for the other species must be calculated using energy adjusted by the lowest energy of each species relative to the lowest energy of vinyl cyanide (peak height for the transition states and well depth for the intermediates).

Beyer–Swinehart Sum and Density of States

Using and from Section 1, we can now compute the sums and densities of states in parallel.

Rate Constants

Using equation (1), we now compute the 13 rate constants, taking care to avoid indeterminate values such as 0/0.

The rate constants can now be visualized. Note that as energy increases, all the reaction rates also increase.

Temporal Dependence of the Various Species

The temporal dependence of the relative concentrations of each stable species (precursor, intermediates and products) can be determined using the rate equations and . The reaction dynamics of the decay of the precursor vinyl cyanide, the subsequent build up and decay of the intermediates, and the formation of products, each normalized with respect to , are plotted at a given internal energy determined by a variable slide bar.

A number of qualitative and quantitative results can be gleaned by simple inspection of this dynamic plot. One example is the amount of time it takes for the reaction to run to completion as energy increases. The plot also indicates that there is never a significant amount of any intermediate apart from that in path III. This is not surprising, given the significantly lower energy for the first transition state in path III relative to the first transition state in pathways IV and V, and given that three pathways share .

A final element determined from this analysis is the exit channel ratios (or product distribution ratios) determined at reaction completion and for a specific system internal energy. This is important because many experimental measurements are only capable of measuring product ratios as a function of precursor internal energy. Thus, experimental scientists require such an RRKM analysis to provide a dynamic and atomistic picture of the chemical reaction. The following code determines and plots the final to ratio at internal energies up to .

This determines the energy at which HNC reaches its maximum relative abundance.

The last plot indicates that production is favored at all system internal energies. Isocyanide () only becomes a noticeable product at a system internal energy above . This is due to paths IV–VII becoming competitive with path I at higher energies. Surprisingly, the product ratio reduces at energies greater than (equivalently, there is a maximum in the product distribution curve). This result would be difficult to predict, even qualitatively, with just the calculated rate constants, demonstrating the usefulness of the powerful tools available in Mathematica.

References

About the Authors

Adam C. Mansell is a Ph.D. candidate in the Bellert Group of the Department of Chemistry and Biochemistry at Baylor University. He holds a B.A. in chemistry from Concordia University Irvine.

David J. Kahle is a computational statistician and assistant professor of statistics at Baylor University. He holds a B.A. in mathematics from the University of Richmond and M.A. and Ph.D. degrees in statistics from Rice University.

Darrin J. Bellert is an associate professor of chemistry at Baylor University. He received his B.S. from Wright State University in Ohio, and his Ph.D. from Florida State University.

Adam C. Mansell
Department of Chemistry and Biochemistry
Baylor University
One Bear Place #97348
Waco, TX 76798
adam_mansell@baylor.edu

David J. Kahle
Department of Statistical Science
Baylor University
One Bear Place #97140
Waco, TX 76798
david.kahle@gmail.com

Darrin J. Bellert
Department of Chemistry and Biochemistry
Baylor University
One Bear Place #97348
Waco, TX 76798
darrin_bellert@baylor.edu