1
Chapter 8. Chemical Dynamics
Chemical dynamics is a field in which scientists study the rates and mechanisms
of chemical reactions. It also involves the study of how energy is transferred among
molecules as they undergo collisions in gas-phase or condensed-phase environments.
Therefore, the experimental and theoretical tools used to probe chemical dynamics must
be capable of monitoring the chemical identity and energy content (i.e., electronic,
vibrational, and rotational state populations) of the reacting species. Moreover, because
the rates of chemical reactions and energy transfer are of utmost importance, these tools
must be capable of doing so on time scales over which these processes, which are often
very fast, take place. Let us begin by examining many of the most commonly employed
theoretical models for simulating and understanding the processes of chemical dynamics.
I. Theoretical Tools for Studying Chemical Change and Dynamics
A. Transition State Theory
The most successful and widely employed theoretical approach for studying reaction
rates involving species that are undergoing reaction at or near thermal-equilibrium
conditions is the transition state theory (TST) of Eyring. This would not be a good way to
model, for example, photochemical reactions in which the reactants do not reach thermal
2
equilibrium before undergoing significant reaction progress. However, for most thermal
reactions, it is remarkably successful.
In this theory, one views the reactants as undergoing collisions that act to keep all of
their degrees of freedom (translational, rotational, vibrational, electronic) in thermal
equilibrium. Among the collection of such reactant molecules, at any instant of time,
some will have enough internal energy to access a transition state (TS) on the Born-
Oppenheimer ground state potential energy surface. Within TST, the rate of progress
from reactants to products is then expressed in terms of the concentration of species that
exist near the TS multiplied by the rate at which these species move through the TS
region of the energy surface.
The concentration of species at the TS is, in turn, written in terms of the equilibrium
constant expression of statistical mechanics discussed in Chapter 7. For example, for a
bimolecular reaction A + B → C passing through a TS denoted AB*, one writes the
concentration (in molecules per unit volume) of AB* species in terms of the
concentrations of A and of B and the respective partition functions as
[AB*] = (qAB*/V)/{(qA/V)( qB/V)} [A] [B].
3
There is, however, one aspect of the partition function of the TS species that is specific to
this theory. The qAB* contains all of the usual translational, rotational, vibrational, and
electronic partition functions that one would write down, as we did in Chapter 7, for a
conventional AB molecule except for one modification. It does not contain a {exp(-hνj
/2kT)/(1- exp(-hνj/kT))} vibrational contribution for motion along the one internal
coordinate corresponding to the reaction path.
Figure 8.1 Typical Potential Energy Surface in Two Dimensions Showing Local Minima,
Transition States and Paths Connecting Them.
In the vicinity of the TS, the reaction path can be identified as that direction along which
the PES has negative curvature; along all other directions, the energy surface is positively
curved. For example, in Fig. 8.1, a reaction path begins at Transition Structure B and is
directed "downhill". More specifically, if one knows the gradients {(?E/?sk) }and
4
Hessian matrix elements { Hj,k = ?2E/?sj?sk}of the energy surface at the TS, one can
express the variation of the potential energy along the 3N Cartesian coordinates {sk} of
the molecule as follows:
E (sk) = E(0) + Σk (?E/?sk) sk + 1/2 Σj,k sj Hj,k sk + …
where E(0) is the energy at the TS, and the {sk} denote displacements away from the TS
geometry. Of course, at the TS, the gradients all vanish because this geometry
corresponds to a stationary point. As we discussed in the Background Material, the
Hessian matrix Hj,k has 6 zero eigenvalues whose eigenvectors correspond to overall
translation and rotation of the molecule. This matrix has 3N-7 positive eigenvalues whose
eigenvectors correspond to the vibrations of the TS species, as well as one negative
eigenvalue. The latter has an eigenvector whose components {sk} along the 3N Cartesian
coordinates describe the direction of the reaction path as it begins its journey from the TS
backward to reactants (when followed in one direction) and onward to products (when
followed in the opposite direction). Once one moves a small amount along the direction
of negative curvature, the reaction path is subsequently followed by taking infinitesimal
“steps” downhill along the gradient vector g whose 3N components are (?E/?sk). Note
that once one has moved downhill away from the TS by taking the initial step along the
negatively curved direction, the gradient no longer vanishes because one is no longer at
the stationary point.
Returning to the TST rate calculation, one therefore is able to express the
concentration [AB*] of species at the TS in terms of the reactant concentrations and a
5
ratio of partition functions. The denominator of this ratio contains the conventional
partition functions of the reactant molecules and can be evaluated as discussed in Chapter
7. However, the numerator contains the partition function of the TS species but with one
vibrational component missing (i.e., qvib = Πk=1,3N-7 {exp(-hνj /2kT)/(1- exp(-hνj/kT))}).
Other than the one missing qvib, the TS's partition function is also evaluated as in Chapter
7. The motion along the reaction path coordinate contributes to the rate expression in
terms of the frequency (i.e., how often) with which reacting flux crosses the TS region
given that the system is in near-thermal equilibrium at temperature T.
To compute the frequency with which trajectories cross the TS and proceed
onward to form products, one imagines the TS as consisting of a narrow region along the
reaction coordinate s; the width of this region we denote δs. We next ask what the
classical weighting factor is for a collision to have momentum ps along the reaction
coordinate. Remembering our discussion of such matters in Chapter 7, we know that the
momentum factor entering into the classical partition function for translation along the
reaction coordinate is (1/h) exp(-ps2/2μkT) dps. Here, μ is the mass factor associated with
the reaction coordinate s. We can express the rate or frequency at which such trajectories
pass through the narrow region of width δs as (ps/μδs), with ps/μ being the speed of
passage (cm s-1) and 1/δs being the inverse of the distance that defines the TS region. So,
(ps/μδs) has units of s-1. In summary, we expect the rate of trajectories moving through
the TS region to be
(1/h) exp(-ps2/2μkT) dps (ps/μδs).
6
However, we still need to integrate this over all values of ps that correspond to enough
energy ps2/2μ to access the TS’s energy, which we denote E*. Moreover, we have to
account for the fact that it may be that not all trajectories with kinetic energy equal to E*
or greater pass on to form product molecules; some trajectories may pass through the TS
but later recross the TS and return to produce reactants. Moreover, it may be that some
trajectories with kinetic energy along the reaction coordinate less than E* can react by
tunneling through the barrier.
The way we account for the facts that a reactive trajectory must have at least E* in
energy along s and that not all trajectories with this energy will react is to integrate over
only values of ps greater than (2μE*)1/2 and to include in the integral a so-called
transmission coefficient κ that specifies the fraction of trajectories crossing the TS that
eventually proceed onward to products. Putting all of these pieces together, we carry out
the integration over ps just described to obtain:
∫ ∫ (1/h) κ exp(-ps2/2μkT) (ps/μδs) ds dps
where the momentum is integrated from ps = (2μE*)1/2 to ∞ and the s-coordinate is
integrated only over the small region δs. If the transmission coefficient is factored out of
the integral (treating it as a multiplicative factor), the integral over ps can be done and
yields the following:
κ (kT/h) exp(-E*/kT).
7
The exponential energy dependence is usually then combined with the partition function
of the TS species that reflect this species’ other 3N-7 vibrational coordinates and
momenta and the reaction rate is then expressed as
Rate = κ (kT/h) [AB*] = κ (kT/h) (qAB* /V)/{(qA/V)( qB/V)} [A] [B].
This implies that the rate coefficient krate for this bimolecular reaction is given in terms of
molecular partition functions by:
krate = κ kT/h (qAB*/V)/{(qA/V)(qB/V)}
which is the fundamental result of TST. Once again we notice that ratios of partition
functions per unit volume can be used to express ratios of species concentrations (in
number of molecules per unit volume), just as appeared in earlier expressions for
equilibrium constants as in Chapter 7.
The above rate expression undergoes only minor modifications when
unimolecular reactions are considered. For example, in the hypothetical reaction A → B
via the TS (A*), one obtains
krate = κ kT/h {(qA*/V)/(qA/V)},
where again qA* is a partition function of A* with one missing vibrational component.
8
Before bringing this discussion of TST to a close, I need to stress that this theory is
not exact. It assumes that the reacting molecules are nearly in thermal equilibrium, so it is
less likely to work for reactions in which the reactant species are prepared in highly non-
equilibrium conditions. Moreover, it ignores tunneling by requiring all reactions to
proceed through the TS geometry. For reactions in which a light atoms (i.e., an H or D
atom) is transferred, tunneling can be significant, so this conventional form of TST can
provide substantial errors in such cases. Nevertheless, TST remains the most widely used
and successful theory of chemical reaction rates and can be extended to include tunneling
and other corrections as we now illustrate.
B. Variational Transition State Theory
Within the TST expression for the rate constant of a bi-molecular reaction, krate = κ
kT/h (qAB*/V)/{(qA/V)(qB/V)}or of a uni-molecular reaction, krate = κ kT/h
{(qA*/V)/(qA/V)}, the height (E*) of the barrier on the potential energy surface
appears in the TS species’ partition function qAB* or qA*, respectively. In particular,
the TS partition function contains a factor of the form exp(-E*/kT) in which the Born-
Oppenheimer electronic energy of the TS relative to that of the reactant species
appears. This energy E* is the value of the potential energy E(S) at the TS geometry,
which we denote S0.
It turns out that the conventional TS approximation to krate over-estimates reaction
rates because it assumes all trajectories that cross the TS proceed onward to products
unless the transmission coefficient is included to correct for this. In the variational
transition state theory (VTST), one does not evaluate the ratio of partition functions
9
appearing in krate at S0, but one first determines at what geometry (S*) the TS partition
function (i.e., qAB* or qA*) is smallest. Because this partition function is a product of
(i) the exp(-E(S)/kT) factor as well as (ii) 3 translational, 3 rotational, and 3N-7
vibrational partition functions (which depend on S), the value of S for which this
product is smallest need not be the conventional TS value S0. What this means is that
the location (S*) along the reaction path at which the free-energy reaches a saddle
point is not the same the location S0 where the Born-Oppenheimer electronic energy
E(S) has its saddle. This interpretation of how S* and S0 differ can be appreciated by
recalling that partition functions are related to the Helmholtz free energy A by q =
exp(-A/kT); so determining the value of S where q reaches a minimum is equivalent
to finding that S where A is at a maximum.
So, in VTST, one adjusts the “dividing surface” (through the location of the
reaction coordinate S) to first find that value S* where krate has a minimum. One then
evaluates both E(S*) and the other components of the TS species partition functions
at this value S*. Finally, one then uses the krate expressions given above, but with S
taken at S*. This is how VTST computes reaction rates in a somewhat different
manner than does the conventional TST. As with TST, the VTST, in the form
outlined above, does not treat tunneling and the fact that not all trajectories crossing
S* proceed to products. These corrections still must be incorporated as an “add-on” to
this theory (i.e., in the κ factor) to achieve high accuracy for reactions involving light
species (recall from the Background Material that tunneling probabilities depend
exponentially on the mass of the tunneling particle).
10
C. Reaction Path Hamiltonian Theory
Let us review what the reaction path is as defined above. It is a path that
i. begins at a transition state (TS) and evolves along the direction of negative curvature on
the potential energy surface (as found by identifying the eigenvector of the Hessian
matrix Hj,k = ?2E/?sk?sj that belongs to the negative eigenvalue);
ii. moves further downhill along the gradient vector g whose components are gk = ?E/?sk’
iii. terminates at the geometry of either the reactants or products (depending on whether
one began moving away from the TS forward or backward along the direction of negative
curvature).
The individual “steps” along the reaction coordinate can be labeled S0, S1, S2, … SP as
they evolve from the TS to the products (labeled SP) and S-R, S-R+1, …S0 as they evolve
from reactants (S-R) to the TS. If these steps are taken in very small (infinitesimal)
lengths, they form a continuous path and a continuous coordinate that we label S.
At any point S along a reaction path, the Born-Oppenheimer potential energy
surface E(S), its gradient components gk(S) = (?E(S)/?sk) and its Hessian components
Hk,j(S) = (?2E(S)/?sk?sj) can be evaluated in terms of derivatives of E with respect to the
3N Cartesian coordinates of the molecule. However, when one carries out reaction path
dynamics, one uses a different set of coordinates for reasons that are similar to those that
arise in the treatment of normal modes of vibration as given in the Background Material.
In particular, one introduces 3N mass-weighted coordinates xj = sj (mj)1/2 that are related
to the 3N Cartesian coordinates sj in the same way as we saw in the Background Material.
11
The gradient and Hessian matrices along these new coordinates {xj} can be
evaluated in terms of the original Cartesian counterparts:
gk’(S) = gk(S) (mk)-1/2
Hj,k’ = Hj,k (mjmk)-1/2.
The eigenvalues {ωk2} and eigenvectors {vk} of the mass-weighted Hessian H’ can then
be determined. Upon doing so, one finds
i. 6 zero eigenvalues whose eigenvectors describe overall rotation and translation of
the molecule;
ii. 3N-7 positive eigenvalues {ωK2} and eigenvectors vK along which the gradient g
has zero (or nearly so) components;
iii. and one eigenvalue ωS2 (that may be positive, zero, or negative) along whose
eigenvector vS the gradient g has its largest component.
The one unique direction along vS gives the direction of evolution of the reaction path (in
these mass-weighted coordinates). All other directions (i.e., within the space spanned by
the 3N-7 other vectors {vK}) possess zero gradient component and positive curvature.
This means that at any point S on the reaction path being discussed
i. one is at a local minimum along all 3N-7 directions {vK} that are transverse to the
reaction path direction (i.e., the gradient direction);
ii. one can move to a neighboring point on the reaction path by moving a small
(infinitesimal) amount along the gradient.
12
In terms of the 3N-6 mass-weighted Hessian’s eigen-mode directions ({vK} and
vS), the potential energy surface can be approximated, in the neighborhood of each such
point on the reaction path S, by expanding it in powers of displacements away from this
point. If these displacements are expressed as components
i. δXk along the 3N-7 eigenvectors vK and
ii. δS along the gradient direction vS,
one can write the Born-Oppenheimer potential energy surface locally as:
E = E(S) + vS δS + 1/2 ωS2 δS2 + ΣK=1,3N-7 1/2 ωK2 δXK2 .
Within this local quadratic approximation, E describes a sum of harmonic potentials
along each of the 3N-7 modes transverse to the reaction path direction. Along the
reaction path, E appears with a non-zero gradient and a curvature that may be positive,
negative, or zero.
The eigenmodes of the local (i.e., in the neighborhood of any point S along the
reaction path) mass-weighted Hessian decompose the 3N-6 internal coordinates into 3N-7
along which E is harmonic and one (S) along which the reaction evolves. In terms of
these same coordinates, the kinetic energy T can also be written and thus the classical
Hamiltonian H = T + V can be constructed. Because the coordinates we use are mass-
weighted, in Cartesian form the kinetic energy T contains no explicit mass factors:
T = 1/2 Σj mj (dsj/dt)2 = 1/2 Σj (dxj/dt)2.
13
This means that the momenta conjugate to each (mass-weighted) coordinate xj, obtained
in the usual way as pj = ?[T-V]/?(dxj/dt) = dxj/dt, all have identical (unit) mass factors
associated with them.
To obtain the working expression for the reaction path Hamiltonian (RPH), one
must transform the above equation for the kinetic energy T by replacing the 3N Cartesian
mass-weighted coordinates {xj} by
i. the 3N-7 eigenmode displacement coordinates δXj,
ii. the reaction path displacement coordinate δS, and
iii. 3 translation and 3 rotational coordinates.
The 3 translational coordinates can be separated and ignored (because center-of-mass
energy is conserved) in further consideration. The 3 rotational coordinates do not enter
into the potential E, but they do appear in T. However, it is most common to ignore their
effects on the dynamics that occurs in the internal-coordinates; this amounts to ignoring
the effects of overall centrifugal forces on the reaction dynamics. We will proceed with
this approximation in mind.
Although it is tedious to perform the coordinate transformation of T outlined
above, it has been done and results in the following form for the RPH:
H = ΣK=1,3N-7 1/2[pK2 + ωK2(S)] + E(S) + 1/2 [pS - ΣK,K’=1,3N-7 pK δXK’ BK,K’]2 /(1+F)
where
14
(1+F) = [1 + ΣK=1,3N-7 δXK BK,S]2.
In the absence of the so-called dynamical coupling factors BK,K’ and BK,S, this expression
for H describes
(1) 3N-7 harmonic-oscillator Hamiltonia 1/2[pK2 + ωK2(S)] each of which has a locally
defined frequency ωK(S) that varies along the reaction path (i.e., is S-dependent);
(2) a Hamiltonian 1/2 pS2 + E(S) for motion along the reaction coordinate S with E(S)
serving as the potential.
In this limit (i.e., with the B factors “turned off”), the reaction dynamics can be simulated
in what is termed a “vibrationally adiabatic” manner by
i. placing each transverse oscillator into a quantum level vK that characterizes the
reactant’s population of this mode;
ii. assigning an initial momentum pS(0) to the reaction coordinate that is
characteristic of the collision to be simulated (e.g., pS(0) could be sampled from a
Maxwell-Boltzmann distribution if a thermal reaction is of interest, or pS(0) could
be chosen equal to the mean collision energy of a beam-collision experiment);
iii. time evolving the S and pS, coordinate and momentum using the above
Hamiltonian, assuming that each transverse mode remains in the quantum state vK
that it had when the reaction began.
The assumption that vK remains fixed, which is why this model is called vibrationally
adiabatic, does not mean that the energy content of the Kth mode remains fixed because
the frequencies ωK(S) vary as one moves along the reaction path. As a result, the kinetic
15
energy along the reaction coordinate 1/2 pS2 will change both because E(S) varies along S
and because ΣK=1,3N-7 hωK2(S) [vK + 1/2] varies along S.
Let’s return now to the RPH theory in which the dynamical couplings among
motion along the reaction path and the modes transverse to it are included. In the full
RPH, the terms BK,K’(S) couple modes K and K’, while BK,S(S) couples the reaction path
to mode K. These couplings are how energy can flow among these various degrees of
freedom. Explicit forms for the BK,K’ and BK,S factors are given in terms of the
eigenvectors {vK, vS} of the mass-weighted Hessian matrix as follows:
BK,K’ = <dvK/dS| vK’>; BK,S = <dvK/dS | vS>
where the derivatives of the eigenvectors {dvK/dS} are usually computed by taking the
eigenvectors at two neighboring points S and S’ along the reaction path:
dvK/dS = {vK(S’) – vK(S))/(S’-S).
In summary, once a reaction path has been mapped out, one can compute, along
this path, the mass-weighted Hessian matrix and the potential E(S). Given these
quantities, all terms in the RPH
H = ΣK=1,3N-7 1/2[pK2 + ωK2(S)] + E(S) + 1/2 [pS - ΣK,K’=1,3N-7 pK δXK’ BK,K’]2 /(1+F)
are in hand.
16
This knowledge can, subsequently, be used to perform the propagation of a set of
classical coordinates and momenta forward in time. For any initial (i.e., t = 0) momenta
pS and pK , one can use the above form for H to propagate the coordinates {δXK, δS} and
momenta {pK, pS} forward in time. In this manner, one can use the RPH theory to follow
the time evolution of a chemical reaction that begins (t = 0) with coordinates and moment
characteristic of reactants under specified laboratory conditions and moves through a TS
and onward to products. Once time has evolved long enough for product geometries to be
realized, one can interrogate the values of 1/2[pK2 + ωK2(S)] to determine how much
energy has been deposited into various product-molecule vibrations and of 1/2 pS2 to see
what the final kinetic energy of the product fragments is. Of course, one also monitors
what fraction of the trajectories, whose initial conditions are chosen to represent some
experimental situation, progress to product geometries vs. returning to reactant
geometries. In this way, one can determine the overall reaction probability.
D. Classical Dynamics Simulation of Rates
One can perform classical dynamics simulations of reactive events without using
the reaction path Hamiltonian. Following a procedure like that outlined in Chapter 7
where condensed-media MD simulations were discussed, one can time-evolve the
Newton equations of motion of the molecular reaction species using, for example, the
Cartesian coordinates of each atom in the system and with either a Born-Oppenheimer
surface or a parameterized functional form. Of course, it is essential for whatever
17
function one uses to accurately describe the reactive surface, especially near the transition
state.
With each such coordinate having an initial velocity (dq/dt)0 and an initial value
q0 , one then uses Newton’s equations written for a time step of duration δt to propagate q
and dq/dt forward in time according, for example , to the following first-order
propagation formula:
q(t+δ t) = q0 + (dq/dt)0 δt
dq/dt (t+δt) = (dq/dt)0 - δt [(?E/?q)0/mq].
Here mq is the mass factor connecting the velocity dq/dt and the momentum pq conjugate
to the coordinate q:
pq = mq dq/dt,
and -(?E/?q)0 is the force along the coordianate q at the “initial” geometry q0. Again, as
in Chapter 7, I should note that the above formulas for propagating q and p forward in
time represent only the most elementary approach to this problem. There are other, more
sophisticated, numerical methods for effecting more accurate and longer-time
propagations, but I will not go into them here. Rather, I wanted to focus on the basics of
how these simulations are carried out.
18
By applying the time-propagation proecess, one generates a set of “new”
coordinates q(t+δt) and new velocities dq/dt(t+δt) appropriate to the system at time t+δt.
Using these new coordinates and momenta as q0 and (dq/dt)0 and evaluating the forces
–(?E/?q)0 at these new coordinates, one can again use the Newton equations to generate
another finite-time-step set of new coordinates and velocities. Through the sequential
application of this process, one generates a sequence of coordinates and velocities that
simulate the system’s dynamical behavior.
In using this kind of classical trajectory approach to study chemical reactions, it is
important to choose the initial coordinates and momenta in a way that is representative of
the experimental conditions that one is attempting to simulate. The tools of statistical
mechanics discussed in Chapter 7 guide us in making these choices. When one attempts,
for example, to simulate the reactive collisions of an A atom with a BC molecule to
produce AB + C, it is not appropriate to consider a single classical (or quantal) collision
between A and BC. Why? Because in any laboratory setting,
1. The A atoms are probably moving toward the BC molecules with a distribution of
relative speeds. That is, within the sample of molecules (which likely contains 1010 or
more molecules), some A + BC pairs have low relative kinetic energies when they
collide, and others have higher relative kinetic energies. There is a probability
distribution P(EKE ) for this relative kinetic energy that must be properly sampled in
choosing the initial conditions.
2. The BC molecules may not all be in the same rotational (J) or vibrational (v) state.
There is a probability distribution function P(J,v) describing the fraction of BC molecules
that are in a particular J state and a particular v state. Initial values of the BC molecule's
19
internal vibrational coordinate and momentum as well as its orientation and rotational
angular momentum must be selected to represent this P(J,v).
3. When the A and BC molecules collide with a relative motion velocity vector v, they do
not all hit "head on". Some collisions have small impact parameter b (the closest distance
from A to the center of mass of BC if the collision were to occur with no attractive or
repulsive forces), and some have large b-values (see Fig. 8.2). The probability function
for these impact parameters is P(b) = 2pi b db, which is simply a statement of the
geometrical fact that larger b-values have more geometrical volume element than smaller
b-values.
Figure 8.2 Coordinates Needed to Characterize an Atom-Diatom Collision Showing the
Impact Parameter b.
20
So, to simulate the entire ensemble of collisions that occur between A atoms and
BC molecules in various J, v states and having various relative kinetic energies EKE and
impact parameters b, one must:
1. run classical trajectories (or quantum propagations) for a large number of J, v, EKE ,
and b values,
2. with each such trajectory assigned an overall weighting (or importance) of
Ptotal = P(EKE ) P(J,v) 2pib db.
After such an ensemble of trajectories representative of an experimental condition
has been carried out, one has available a great deal of data. This data includes knowledge
of what fraction of the trajectories produced final geometries characteristic of products,
so the net reaction probability can be calculated. In addition, the kinetic and potential
energy content of the internal (vibrational and rotational) modes of the product molecules
can be interrogated and used to compute probabilities for observing products in these
states. This is how classical dynamics simulations allow us to study chemical reactions
and/or energy transfer.
E. RRKM Theory
Another theory that is particularly suited for studying uni-molecular decomposition
reactions is named after the four scientists who developed it- Rice, Ramsperger, Kassel,
and Marcus. To use this theory, one imagines an ensemble of molecules that have been
21
“activated” to a state in which they possess a specified total amount of internal energy E
of which an amount E*rot exists as rotational energy and the remainder as internal
vibrational energy.
The mechanism by which the molecules become activated could involve collisions
or photochemistry. It does not matter as long as enough time has passed to permit one to
reasonably assume that these molecules have the energy E-E*rot distributed randomly
among all their internal vibrational degrees of freedom. When considering thermally
activated unimolecular decomposition of a molecule, the implications of such
assumptions are reasonably clear. For photochemically activated unimolecular
decomposition processes, one usually also assumes that the molecule has undergone
radiationless relaxation and returned to its ground electronic state but in a quite
vibrationally “hot” situation. That is, in this case, the molecule contains excess
vibrational energy equal to the energy of the optical photon used to excite it. Finally,
when applied to bimolecular reactions, one assumes that collision between the two
fragments results in a long-lived complex. The lifetime of this intermediate must be long
enough to allow the energy E-E*rot, which is related to the fragments’ collision energy, to
be randomly distributed among all vibrational modes of the collision complex.
For bimolecular reactions that proceed “directly” (i.e., without forming a long-lived
intermediate), one does not employ RRKM-type theories because their primary
assumption of energy randomization almost certainly would not be valid in such cases.
The RRKM expression of the unimolecular rate constant for activated molecules A*
(i.e., either a long-lived complex formed in a bimolecular collision or a “hot” molecule)
dissociating to products through a transition state, A* → TS → P, is
22
krate = G(E-E0 –E’tot )/(N(E-E*rot) h).
Here, the total energy E is related to the energies of the activated molecules by
E = E*rot + E*vib
where E*rot is the rotational energy of the activated molecule and E*vib is the vibrational
energy of this molecule. This same energy E must, of course, appear in the transition state
where it is decomposed as an amount E0 needed to move from A* to the TS (i.e, the
energy needed to reach the barrier) and vibrational (E'vib) , translational (E'trans along the
reaction coordiate), and rotational (E'rot) energies:
E = E0 + E’vib + E’trans + E’rot .
In the rate coefficient expression, G(E-E0 –E’rot ) is the total sum of internal
vibrational quantum states that the transition state possesses having energies up to and
including E-E0 –E’rot . This energy is the total energy E but with the activation energy E0
removed and the overall rotational energy E’rot of the TS removed. The quantity
N(E-E*rot) is the density of internal vibrational quantum states (excluding the mode
describing the reaction coordinate) that the activated molecule possesses having an
energy between E-E*rot and E-E*rot + dE. In this expression, the energy E-E*rot is the total
energy E with the rotational energy E*rot of the activated species removed.
23
In the most commonly employed version of RRKM theory, the rotational energies of
the activated molecules E*rot and of the TS E’rot are assumed to be related by
E*rot – E’rot = J(J+1) h2/8pi2 {1/I* - 1/I’} = E*rot {1 – I*/I’}.
Here I* and I’ are the average (taken over the three eigenvalues of the moment inertia
tensors) moments of inertia of the activated molecules and TS species, respectively. The
primary assumption embodied in the above relationship is that the rotational angular
momenta of the activated and TS species are the same, so their rotational energies can be
related, as expressed in the equation, to changes in geometries as reflected in their
mometns of inertia. Because RRKM theory assumes that the vibrational energy is
randomly distributed, its fundamental rate coefficient equation
krate = G(E-E0 –E’tot )/(N(E-E*rot) h) depends on the total energy E, the energy E0 required
to access the TS, and the amount of energy contained in the rotational degrees of freedom
that is thus not available to the vibrations.
To implement a RRKM rate coefficient calculation, one must know
(i) the total energy E available,
(ii) the barrier energy E0,
(iii) the geometries (and hence the moments of inertia I* and I’) of the activated
molecules and of the TS, respectively,
(iv) the rotational energy E*rot of the activated molecules, as well as
(v) all 3N-6 vibrational energies of the activated molecules and all 3N-7 vibrational
energiers of the TS (i.e., excluding the reaction coordinate).
24
The rotational energy of the TS species can then be related to that of the activated
molecules through E*rot – E’rot = E*rot {1 – I*/I’}.
To simulate an experiment in which the activated molecules have a thermal
distribution of rotational energies, the RRKM rate constant is computed for a range of
E*rot values and then averaged over E*rot using the thermal Boltzmann population
(2J+1) exp(-J(J+1)h2/(8pi2I*kT)) as a weighting factor. This can be carried out, for
example, using the MC process for selecting rotational J values. This then produces a rate
constant for any specified total energy E. Alternatively, to simulate experiments in which
the activated species are formed in bimolecular collisions at a specified energy E, the
RRKM rate coefficient is computed for a range of E*rot values with each E*rot related to
the collisional impact parameter b that we discussed earlier. In that case, the collisional
angular momentum J is given as J = μ v b, where v is the collision speed (related to the
collision energy) and μ is the reduced mass of the two colliding fragments. Again using
E*rot – E’rot = E*rot {1 – I*/I’} the TS rotational energy can be related to that of the
activated species. Finally, the RRKM rate coefficient is evaluated by averaging the result
over a series of impact parameters b (each of which implies a J value and thus an E*rot )
with 2pib db as the weighting factor.
The evaluation of the sum of states G(E-E0 –E’tot ) and the density of states
N(E-E*rot) that appear in the RRKM expression is usually carried out using a state-
counting algorithm such as that implemented by Beyer and Swinehart (Commun. Assoc.
Comput. Machin. 16, 372 (1973)). This algorithm uses knowledge of the 3N-6 harmonic
vibrational frequencies of the activated molecules and the 3N-7 frequencies of the TS and
determines how many ways a given amount of energy can be distributed among these
25
modes. By summing over all such distributions for energy varying from zero to E, the
algorithm determines G(E). By taking the difference G(E+δE) – G(E), it determines
N(E)δE.
F. Correlation Function Expressions for Rates
Recall from Chapter 6 that rates of photon absorption can, in certain
circumstances, be expressed either (a) in terms of squares of transition dipole matrix
elements connecting each initial state Φi to each final state Φf,
| E0 ? <Φf | μ | Φi> |2
or (b) in terms of the equilibrium average of the product of a transition dipole vector at
time t=0 dotted into this same vector at another time t
Σi ρi <Φi | E0 ? μ E0 ? μ (t) | Φi>
That is, these rates can be expressed either in a state-to-state manner or in a time-
dependent correlation function framework. In Chapter 7, this same correlation function
approach was examined further.
In an analogous fashion, it is possible to express chemical reaction rate constants
in a time-domain language again using time correlation functions. The TST (or VTST)
and RRKM expressions for the rate constant krate all involve, through the partition
26
functions or state densities, the reactant and transition-state energy levels and
degeneracies. These theories are therefore analogs of the state-to-state photon-absorption
rate equations.
To make the connection between the state-to-state and time-correlation function
expressions, one can begin with a classical expression for the rate constant given below:
Here Qr is the partition function of the reactant species, L is the number of coordinates
and momenta upon which the Hamiltonian H(p,q) depends, and β is 1/kT. The flux factor
F and the reaction probability χ are defined in terms of a dividing surface which could,
for example, be a plane perpendicular to the reaction coordinate S and located along
the reaction path that was discussed earlier in this Chapter in Section IC. Points on such
a surface can be defined by specifying one condition that the L coordinates {qj} must
obey, and we write this condition as
f(q) = 0.
Points lying where f(q) < 0 are classified as lying in the reactant region of coordinate
space, while those lying where f > 0 are in the product region. For example, if the
dividing surface is defined as mentioned above as being a plane perpendicular to the
reaction path, the function f can be written as:
k(T) = Q
r
?1
(2pih)
? L
dpdqe
? βH(q, p)
∫ F(p,q)χ(p,q)
27
f(q) = (S - S0).
Here, S is the reaction coordinate (which, of course, depends on all of the q variables)
and S0 is the value of S at the dividing surface. If the dividing surface is placed at the
transition state on the energy surface, S0 vanishes because the transition state is then, by
convention, the origin of the reaction coordinate.
So, now we see how the dividing surface can be defined, but how are the flux
F and probability χ constructed? The flux factor F is defined in terms of the dividing
surface function f(q) as follows:
F(p,q) = d h(f(q))/dt
= (dh/df) (df/dt)
=(dh/df) Σj ?f/?qj (dqj/dt)
= δ(f(q)) Σj ?f/?qj (dqj/dt).
Here, h(f(q)) is the Heaviside step function (h(x) = 1 if x>0; h(x) = 0 if x < 0), whose
derivative dh(x)/dx is the Dirac delta function δ(x), and the other identities follow by
using the chain rule. When the dividing surface is defined in terms of the reaction path
coordinate S as introduced earlier (i.e., f(q) = (S - S0)), the factor Σj ?f/?qj (dqj/dt)
contains only one term when the L coordinates {qj} are chosen, as in the reaction path
28
theory, to be the reaction coordinate S and L-1 coordinates {q’j} = q’ perpendicular to the
reaction path. For such a choice, one obtains
Σj ?f/?qj (dqj/dt) = dS/dt = PS/mS
where PS is the momentum along S and mS is the mass factor associated with S in the
reaction path Hamiltonian. So, in this case, the total flux factor F reduces to:
F(p,q) = δ(S-S0) PS/mS.
We have seen exactly this construct before in Section IA where the TST expression for
the rate coefficient was developed.
The reaction probability factor χ(p,q) is defined in terms of those trajectories that
evolve, at long time t → ∞, onto the product side of the dividing surface; such trajectories
obey
χ(p,q) = limt → ∞ h(f(q(t))) = 1.
This long-time limit can, in turn, be expressed in a form where the flux factor again
occurs:
lim
t →∞
h( f (q(t))) =
dh( f (q(t)))
dt
0
∞
∫ dt = Fdt
0
∞
∫
29
In this expression, the flux F(t) pertains to coordinates q(t) and momenta p(t) at t > 0.
Because of time reversibility, the integral can be extended to range from from t = - ∞ to t
= ∞.
Using the expressions for χ and for F as developed above in the equation for the
rate coefficient given at the beginning of this Section allows the rate coefficient k(T) to
be rewritten as follows:
In this form, the rate constant k(T) appears as an equilibrium average (represented by the
integral over the initial values of the variables p and q with the Qr-1 (2pih)-L exp(-βH)
weighting factor) of the time correlation function of the flux F:
To evaluate the rate constant in this time-domain framework for a specific
chemical reaction, one would proceed as follows.
i. Run an ensemble of trajectories whose initial coordinates and momenta {q.p} are
selected (e.g., using Monte-Carlo methods discussed in Chapter 7) from a distribution
with exp(-βH) as its weighting factor.
k(T) = Q
r
?1
(2pih)
? L
dpdqe
? βH(q, p)
∫ F(p,q)χ(p,q)
= Q
r
?1
(2pih)
? L
dt
?∞
∞
∫ dpdqe
?βH (q, p)
∫ F( p,q)F( p(t),q(t))
dt
?∞
∞
∫ F(p,q)F( p(t),q(t))
30
ii. Make sure that the initial coordinates {q} lie on the dividing surface because the flux
expression contains the δ(f(q)) factor;
iii. Monitor each trajectory to observe when it again crosses the dividing surface (i.e.,
when {q(t)} again obeys f(q(t)) = 0; at which time the quantity
iv. F(p(t),q(t)) can be evaluated as F(p,q) = δ(f(q)) Σj ?f/?qj (dqj/dt), using the coordinates
and momenta at time t to compute these quantities.
Using a planar dividing surface attached to the reaction path at S = S0 as noted
earlier allows F(q,p) to be calculated in terms of the initial (t=0) momentum lying along
the reaction path direction as , F(p,q) = δ(S-S0) PS/mS and permits F(p(t),q(t)) to be
computed when the trajectory again crosses this surface at at time t as F(p(t),q(t)) = δ(S-
S0) PS(t)/mS. So, all that is really needed if the dividing surface is defined in this manner
is to start trajectories with S = S0; to keep track of the initial momentum along S; to
determine at what times t the trajectory returns to S = S0; and to form the product (PS/mS)
(PS(t)/mS) for each such time. It is in this manner that one can compute flux-flux
correlation functions and, thus, the rate coefficient.
Notice that trajectories that undergo surface re-crossings contribute negative
terms to the flux-flux correlation function computed as discussed above. That is, a
trajectory with a positive initial value of (PS/mS) can, at some later time t, cross the
dividing surface with a negative value of (PS(t)/mS) (i.e, be directed back toward
reactants). This re-crossing will contribute a negative value, via. the product (PS/mS)
(PS(t)/mS), to the total correlation function, which integrates over all times. Of course, if
this same trajectory later undergoes yet another crossing of the dividing surface at t' with
31
positive PS(t'), it will contribute a positive term to the correlation function via. (PS/mS)
(PS(t')/mS). Thus, the correlation function approach to computing the rate coefficient can
properly account for surface re-crossings, unlike the TST which requires one to account
for such effects (and tunneling) in the transmission coefficient κ.
G. Wave Packet Propagation
The preceding discussion should have made it clear that it is very difficult to
propagate wave functions rigorously using quantum mechanics. On the other hand, to
propagate a classical trajectory is relatively straightforward. There exists a powerful tool
that allows one to retain much of the computational ease and convenient interpretation of
the classical trajectory approach while also incorporating quantum effects that are
appropriate under certain circumstances. In this wave packet propagation approach, one
begins with a quantum mechanical wave function that is characterized by two parameters
that give the average value of the position and momentum along each coordinate. One
then propagates not the quantum wave function but the values of these two parameters,
which one assumes will evolve according to Newtonian dynamics. Let's see how these
steps are taken in more detail and try to understand when such an approach is expected to
work or to fail.
First, the form of the so-called wave packet quantum function is written as
follows:
Ψ(q,Q, P) = ΠJ=1,N (2pi<δqJ2>)-1/2 exp[(iPJ qJ/ h) -(qJ - QJ)2 /4<δqJ2>].
32
Here, we have a total of N coordinates that we denote {qJ : J = 1, N}. It is these
coordinates that the quantum wave function depends upon. The total wave function is a
product of terms, one for each coordinate. Notice that this wave function has two distinct
ways in which the coordinate qJ appear. First, it has a Gaussian (exp[- (qJ - QJ)2 /4<δqJ2>])
dependence centered at the values QJ and having Gaussian width factors related to <qJ2>.
This dependence tends to make the wave function's amplitude largest when qJ is close to
QJ. Secondly, it has a form exp[(iPJ qJ/ h)] that looks like the travelling wave that we
encountered in the Background Material in which the coordinate qJ moves with
momentum PJ . So, these wave packet functions have built into them characteristics that
allow them to describe motion (via. the PJ) of an amplitude that is centered at QJ with a
width given by the parameter <qJ2>.
The parameters PJ and QJ we assume in this approach to chemical dynamics will
undergo classical time evolution according to the Newton equations:
dQJ/dt = PJ/mJ
dPJ/dt = - ?E/?QJ
where E is the potential energy surface (Born-Oppenheimer or force field) upon which
we wish to propagate the wave packet, and mJ is the mass associated with coordinate qJ .
The QJ and PJ parameters can be shown to be the expectation values of the coordinates qJ
and momenta -i h?/?qJ for the above function:
33
QJ = ∫ Ψ* qJ Ψ dq,
PJ = ∫ Ψ* (- i h ?/?qJ) Ψ dq.
Moreover, the <qJ2> parameter appearing in the Gaussian part of the function can be
shown to equal the dispersion or "spread" of this wave function along the coordinate qJ:
<qJ2> = ∫ Ψ* (qJ - QJ)2 Ψ dq.
There is an important characteristic of the above Gaussian wave packet functions
that we need to point out. It turns out that functions of this form:
Ψ(q,Q(t), P(t)) = ΠJ=1,N (2pi<δqJ2>)-1/2 exp[(iPJ(t) qJ/ h) -(qJ - QJ(t))2 /4<δqJ2>]
can be shown to have uncertainties in qJ and in - i h ?/?qJ that are as small as possible:
<(qJ –QJ)2> <(- i h ?/?qJ – PJ)2> = h2/4.
The Heisenberg uncertainty relation, which is discussed in many texts dealing with
quantum mechanics, says that this product of coordinate and momentum dispersions must
be greater than or equal to h2/4. In a sense, the Gaussian wave packet function is the most
34
classical function that one can have because its uncertainty product is as small as possible
(i.e, equals h2/4). We say this is the most classical possible quantum function because in
classical mechanics, both the coordinate and the momentum can be known precisely. So,
whatever quantum wave function allows these two variables to be least uncertain is the
most classical.
To use wave packet propagation to simulate a chemical dynamics event, one
begins with a set of initial classical coordinates and momenta {QJ(0), PJ(0)} as well as a
width <qJ2> or uncertainty for each coordinate. Each width must be chosen to represent
the range of that coordinate in the experiment that is to be simulated. For example,
assume one were to represent the dynamics of a wave function that is prepared by photon
absorption of a v = 0 vibrational state of the H-Cl molecule from the ground 1Σ state to an
excited-state energy surface (E(R)). Such a situation is described qualitatively in Fig. 8.3.
In this case, one could choose <δR2> to be the half width of the v = 0 harmonic (or
Morse) oscillator wave function χ0(R) of H-Cl, and take P(0) = 0 (because this is the
average value of the momentum for χ0) and R(0) = Req, the equilibrium bond length.
For such initial conditions, classical Newtonian dynamics would then be used to
propagate the QJ and PJ. In the H-Cl example, introduced above, this propagation would
be performed using the excited-state energy surface for E since, for t > 0, the molecule is
assumed to be on this surface. The total energy at which the initial wave packet it
delivered to the upper surface would be dictated by the energy of the photon used to
perform the excitation. In Fig. 8.3, two such examples are shown.
Once the packet is on the upper surface, its position Q and momentum P begin to
change according to the Newton equations. This, in turn, causes the packet to move as
35
shown for several equally spaced time “steps” in Fig. 8.3 for the two different photons’
cases. At such subsequent times, the quantum wave function is then assumed, within this
model, to be given by:
Ψ(q,Q(t), P(t)) = ΠJ=1,N (2pi<δqJ2>)-1/2 exp[(iPJ(t) qJ/ h) -(qJ - QJ(t))2 /4<δqJ2>].
That is, it is taken to be of the same form as the initial wave function but to have simply
moved its center from Q(0) to Q(t) with a momentum that has changed from P(0) to P(t).
36
Ground-State Potential
Excited-State Potential
Initital Wave Packet
Wave Packet Just After Absorption of Low-Energy
Photon
Wave Packet Just After Absorption of Higher-Energy
Photon
t=0 t1 t2
t3
t=0 t1 t2
t3 t4 t5
Figure 8.3 Propagation of Wave Packet Prepared by Absorption of Two Different
Photons.
It should be noticed that the time evolution of the wave packet shown in Fig. 8.3 displays
clear classical behavior. For example, as time evolves, it moves to large R-values and its
speed (as evidenced by the spacings between neighboring packets for equal time steps) is
37
large when the potential is low and small when the potential is higher. As we learned in
the Background Material and in Chapter 6, the time correlation function
C(t) = <Ψ(q,Q(0),P(0)) | Ψ(q,Q(t),P(t))>
can be used to extract spectral information by Fourier transformation. For the H-Cl
example considered here, this correlation function will be large at t = 0 but will decay in
magnitude as the wave packet Ψ(q,Q(t),P(t)) moves to the right (at t1, t2, etc.) because its
overlap with Ψ(q,Q(0),P(0)) becomes smaller and smaller as time evolves. This decay in
C(t) will occur more rapidly for the high-energy photon case because Ψ(q,Q(t),P(t))
moves to the right more quickly because the classical momentum P(t) grows more
rapidly. These dynamics will induce exponential decays in C(t) (i.e., C(t) will vary as
exp(-t/τ1)) at short times.
In fact, the decay of C(t) discussed above produces, when C(t) is Fourier
transformed, the primary characteristic of the correlation function for the higher-energy
photon case where dissociation ultimately occurs. In such photo-dissociation spectra, one
observes a Lorentzian line shape whose width is characterized by the decay rate (1/τ1),
which, in turn, relates to the total energy of the packet and the steepness of the excited-
state surface. This steepness determines how fast P(t) grows, which then determines how
fast the H-Cl bond fragments.
In the lower-energy photon case shown in Fig. 8.3, a qualitatively different
behavior occurs in C(t) and thus in the spectrum. The packet’s movement to larger R
causes C(t) to initially undergo exp(-t/τ1) decay. However, as the packet moves to its
38
large-R turning point (shortly after time t3), it strikes the outer wall of the surface where
it is reflected. Subsequently, it undergoes motion to smaller R, eventually returning to its
initial location. Such recurrences, which occur on time scales that we denote τ2, are
characteristic of bound motion in contrast to the directly dissociative motion discussed
earlier. This recurrence will cause C(t) to again achieve a large amplitude, but, C(t) will
subsequently again undergo exp(-t/τ1) decay as the packet once again departs. Clearly, the
correlation function will display a series of recurrences followed by exponential decays.
The frequency of the recurrences is determined by the frequency with which the packet
traverses from its inner to outer turning points and back again, which is proportional to
1/τ2. This, of course, is the vibrational period of the H-Cl bond. So, in such bound-motion
cases, the spectrum (i.e., the Fourier transform of C(t)) will display a series of peaks
spaced by (1/τ2) with the envelope of such peaks having a width determined by 1/τ1.
In more complicated multi-mode cases (e.g., in molecules containing several
coordinates), the periodic motion of the wave packet usually shows another feature that
we have not yet discussed. Let us, for simplicity, consider a case in which only two
coordinates are involved. For the wave packet to return to (or near) its initial location
enough time must pass for both coordinates to have undergone an excursion to their
turning points and back. For example, consider the situation in which one coordinate’s
vibrational frequency is ca. 1000 cm-1 and the other’s is 300 cm-1; these two modes then
require ca. 1/30 ps and 1/9 ps, respectively, to undergo one complete oscillation. At t = 0,
the wave packet, which is a product of two packets, ΠJ=1,2 (2pi<δqJ2>)-1/2 exp[(iPJ(t) qJ/ h) -
(qJ - QJ(t))2 /4<δqJ2>], one for each mode, produces a large C(t). After 1/30 ps, the first
mode’s coordinate has returned to its initial location, but the second mode is only 9/30 of
39
the way along in its periodic motion. Moreover, after 1/9 ps, the second mode’s
coordinate has returned to near where it began, but now the first mode has moved away.
So, at both 1/30 ps and 1/9 ps, the correlation function will not be large because one of
the mode contribution to C(t) = <Ψ(q,Q(0),P(0)) | Ψ(q,Q(t),P(t))> will be small.
However, after 1/3 ps, both modes’ coordinates will be in positions to produce a large
value of C(t); the high-frequency mode will have undergone 10 oscillations, and the
lower-frequency mode will have undergone 3 oscillations. My point in discussing this
example has been to illustrate that molecules having many coordinates can produce
spectra that display rather complicated patterns but which, in principle, can be related to
the time evolution of these coordinates using the correlation function’s connection to the
spectrum.
Of course, there are problems that arise in using the wave packet function to
describe the time evolution of a molecule (or any system that should be treated using
quantum mechanics). One of the most important limitations of the wave packet approach
to be aware of relates to it inability to properly treat wave reflections. It is well know that
when a wave strikes a hard wall it is reflected by the wall. However, when, for example,
a water wave moves suddenly from a region of deep water to a much more shallow
region, one observes both a reflected and a transmitted wave. In the discussion of
tunneling resonances given in the Background Material, we also encountered reflected
and transmitted waves. Furthermore, when a wave strikes a barrier that has two or more
holes or openings in it, one observes wave fronts coming out of these openings. The
problem with the most elementary form of wave packets presented above is that each
packet contains only one “piece”. It therefore can not break into two or more “pieces” as
40
it, for example, reflects from turning points or passes through barriers with holes.
Because such wave packets can not fragment into two or more packets that subsequently
undergo independent dynamical evolution, they are not able to describe dynamical
processes that require multiple-fragmentation events. It is primarily for this reason that
wave packet approaches to simulating dynamics are usually restricted to treating short-
time dynamics where such fragmentation of the wave packet is less likely to occur.
Prompt molecular photo-dissociation processes such as we discussed above is a good
example of such a short-time phenomenon.
H. Surface Hopping Dynamics
There are, of course, chemical reactions and energy transfer collisions in which
two or more Born-Oppenheimer energy surfaces are involved. Under such circumstances,
it is essential to have available the tools needed to describe the coupled electronic and
nuclear-motion dynamics appropriate to this situation.
The way this problem is addressed is by returning to the Schr?dinger equation
before the single-surface Born-Oppenheimer approximation was made and expressing the
electronic wave function Ψ(r|R), which depends on the electronic coordinates {r} and
the nuclear coordinates {R}, as:
Ψ(r|R) = ΣJ aJ(t) ψJ(r|R) exp[-i/ h ∫ HJ,J (R) dt].
41
Here, ψJ(r|R) is the Born-Oppenheimer electronic wave function belonging to the Jth
electronic state, and HJ,J(R) is the expectation value of H for this state. The aJ(t) are
amplitudes that will eventually relate to the probability that the system is “on” the Jth
surface. Substituting this expansion into the time-dependent Schr?dinger equation
i h ?Ψ/?t = H Ψ
followed by multiplying on the left by ψ*K(r|R) and integrating over the electronic
coordinates {r}, gives an equation for the aK(t) amplitudes:
i h daK/dt = ΣJ≠K { HK,J - i h <ψK|dψJ/dt>} exp[-i/ h ∫ (HJ,J (R) – HK,K(R)) dt].
Here, HK,J is the Hamiltonian matrix that couples ψK to ψJ. The integral appearing in the
exponential is taken from an initial time ti when one is assumed to know that the system
resides on a particular Born-Oppenheimer surface (K) up to the time t, at which the
amplitude for remaining on this surface aK as well as the amplitudes for hopping to other
surfaces {aJ} are needed. This differential equation for the amplitudes is solved
numerically by starting at ti with aK = 1 and aJ≠K = 0 and propagating the amplitudes’
values forward in time.
The next step is to express <ψK|dψJ/dt>, using the chain rule, in terms of
derivatives with respect to the nuclear coordinates {R} and the time rate of change of
these coordinates:
42
<ψK|dψJ/dt> = Σa <ψK|dψJ/dRa> dRa/dt.
This is how a classical dynamical treatment of the nuclear motions will be introduced
(i.e., by assuming that the nuclear coordinates Ra and dRa/dt obey Newton equations). So,
now the equations for the aK(t) read as follows:
i h daK/dt = ΣJ≠K { HK,J - i h Σa dRa/dt <ψK|dψJ/dRa>} exp[-i/ h ∫ (HJ,J (R) – HK,K(R)) dt].
The <ψK|dψJ/dRa> are called non-adiabatic coupling matrix elements, and it is their
magnitudes that play a central role in determining how efficient surface hoppings are.
These matrix elements are becoming more commonly available in widely utilized
quantum chemistry and dynamics computer packages (although their efficient evaluation
remains a challenge that is undergoing significant study).
In addition to the above prescription for calculating amplitudes (the probabilities
then being computed as |aJ|2), one also needs to identify (using, perhaps the kind of
strategy discussed in Chapter 3) the seam at which the surfaces of interest intersect. Let
us, for example, assume that there are two surfaces that undergo such an intersection in a
subspace of the 3N-6 dimensional energy surfaces H1,1 and H2,2, and let us denote nuclear
coordinates {Ra} that lie on this seam as obeying
F(Ra) = 0.
following the discussion of Chapter 3.
43
To utilize the most basic form of surface hopping theory, one proceeds as follows:
1. One begins with initial values of the nuclear coordinates {Ra} and their velocities
{dRa/dt} that properly characterize the kind of collision or reaction one wishes to
simulate. Of course, one may have to perform many such surface hopping trajectories
with an ensemble of initial conditions chosen to properly describe such an experimental
situation. In addition, one specifies which electronic surface (say the Kth surface) the
system is initially on.
2. For each such set of initial conditions, one propagates a classical trajectory describing
the time evolution of the {Ra} and {dRa/dt} on this initial surface, until such a trajectory
approaches a geometry lying on the intersection seam.
3. As one is propagating the classical trajectory from ti up to the time t1 when it
approaches the seam, one also propagates the (in this example, two) differential equations
i h daK/dt = ΣJ≠K { HK,J - i h Σa dRa/dt <ψK|dψJ/dRa>} exp[-i/ h ∫ (HJ,J (R) – HK,K(R)) dt]
that produce the time evolution of the {aJ} amplitudes.
4. Once the classical trajectory reaches the seam, the initial trajectory is halted and two
(in the case of M coupled surfaces, M) new trajectories are initiated, one beginning (at t1)
on each of the (two in our case) coupled surfaces. Each of these new trajectories (two )
are assigned a probability weighting factor computed as the square of the amplitude
belonging to the surface upon which the trajectory now resides: |aJ|2.
5. Each of the new trajectories is subsequently allowed to propagate from t1 (where the
nuclear coordinates are what they were just prior to diving the one initial trajectory into
two (or M) trajectories) but now on different Born-Oppenheimer surfaces. This
44
propagation is continued until one of these trajectories again approaches a geometry on
the intersection seam, at which time the trajectory hopping scheme is repeated.
Clearly, it is possible that the one initial trajectory will generate a large number of
subsequent trajectories each of which is assigned an overall weighting factor taken to be
the product of all of the |aJ|2 probabilities associated with trajectories from which this
particular trajectory was descended. It is these overall weightings that one uses to predict,
one this process has been carried out long enough to allow the reaction or collision to
proceed to completion, the final yield of products in each Born-Oppenheimer state. This
surface-hopping algorithm remains one of the most widely used approaches to treating
such coupled-state dynamics.
II. Experimental Probes of Reaction Dynamics
A. Spectroscopic Methods
To follow the rate of any chemical reaction, one must have a means of monitoring
the concentrations of reactant or product molecules as time evolves. In the majority of
current experiments that relate to reaction dynamics, one uses some form of
spectroscopic or alternative physical probe (e.g., an electrochemical signature) to monitor
45
these concentrations as functions of time. Of course, in all such measurements, one must
know how the intensity of the signal detected relates to the concentration of the
molecules that cause the signal. For example, in most absorption experiments, as
illustrated in Fig. 8.4, light is passed through a sample of “thickness” L and the intensity
of the light beam in the absence of the sample I0 and with the sample present I are
measured.
Figure 8.4 Typical Beer’s –Law Experiment in Which a Light Beam of Intensity I0 is
Passed Through a Sample of Thickness L.
The Beer-Lambert law:
log(I0/I) = ε [A] L
Intensity ofIncident Light I
0
Intensity ofTransmitted Light I
Sample of thickness L
46
then allows the concentration [A] of the absorbing molecules to be determined, given the
path length L over which absorption occurs and given the extinction coefficient ε of the
absorbing molecules.
These extinction coefficients, which relate to the electric dipole matrix elements
that are discussed in Chapter 6, are usually determined empirically by preparing a known
concentration of the absorbing molecules and measuring the I0/I ratio that this
concentration produces in a cell of length L. For molecules and ions that are extremely
reactive, this “calibration” approach to determining ε is often not feasible because one
can not prepare a sample with a known and stable concentration. In such cases, one often
must resort to using the theoretical expressions given in Chapter 6 (and discussed in most
textbooks on molecular spectroscopy) to compute ε in terms of the wave functions of the
absorbing species. In any event, one must know how the strength of the signal relates to
the concentrations of the species if one wishes to monitor chemical reaction or energy
transfer rates.
Because modern experimental techniques are capable of detecting molecules in
particular electronic and vibration-rotation states, it has become common to use such
tools to examine chemical reaction dynamics on a state-to-state level and to follow
energy transfer processes, which clearly require such state-specific data. In such
experiments, one seeks to learn the rate at which reactants in a specific state Φi react to
produce products in some specific state Φf. One of the most common ways to monitor
such state-specific rates is through a so-called “pump-probe” experiment in which
i. A short-duration light source is used to excite reactant molecules to some
specified initial state Φi. Usually a tunable laser is used because its narrow
47
frequency spread allows specific states to be pumped. The time at which this
“pump” laser thus prepares the excited reactant molecules in state Φι defines t = 0.
ii. After a “delay time” of duration τ, a second light source is used to “probe” the
product molecules that have been formed in various final states, Φf. Usually, the
frequency of this probe source is scanned so that one can examine populations of
many such final states.
The concentrations of reactant and products molecules in the initial and final states, Φi
and Φf, are determined by the Beer-Lambert relation assuming that the extinction
coefficients εi and εf for these species and states absorption are known. In the former case,
the extinction coefficient ει relates to absorption of the pump photons to prepare reactant
molecules in the specified initial state. In the latter, εf refers to absorption of the product
molecules that are created in the state Φf. Carrying out a series of such final-state
absorption measurements at various delay times τ allows one to determine the
concentration of these states as a function of time.
This kind of laser pump-probe experiment is used not only to probe specific
electronic or vibration/rotation states of the reactants and products but also when the
reaction is fast (i.e., complete in 10-4s or less). In these cases, one is not using the high
frequency resolution of the laser but its fast time response. Because laser pulses of quite
short duration can be generated, these tools are well suited in such fast chemical reaction
studies. The reactions can be in the gas phase (e.g., fast radical reactions in the
atmosphere or in explosions) or in solution (e.g., photo-induced electron transfer
reactions in biological systems).
48
B. Beam Methods
Another approach to probing chemical reaction dynamics is to use a beam of reactant
molecules A that collides with other reactants B that may also in a beam or in a “bulb” in
equilibrium at some temperature T. Such crossed-beam and beam-bulb experiments are
illustrated in Fig. 8.5.
49
Figure 8.5 Typical Crossed-Beam and Beam-Bulb Experimental Setups.
Almost always, these beam and bulb samples contain molecules, radicals, or ions in the
gas phase, so these techniques are most prevalent in gas-phase dynamics studies.
The advantages of the crossed-beam type experiments are that:
(i) one can control the velocities, and hence the collision energies, of both reagents,
Beam of A molecules
moving with velocity
v.
Product Cmolecules
detected exitingwith velocity v'
in internal state jat angle θθ
Bulb containingB molecules at
temperature T
Beam of A molecules
moving with velocity
v.
Product Cmolecules
detected exitingwith velocity v''
in internal state jat angle θθ
Beam of B molecules
moving with velocity
v'.
Collision chamberin which A and B
collide and react
50
(ii) one can examine the product yield as a function of the angle θ through which the
products are scattered,
(iii) one can probe the velocity of the products and,
(iv) by using spectroscopic methods, one can determine the fraction of products
generated in various internal (electronic/vibrational/rotational) states.
Such measurements allow one to gain very detailed information about how the
reaction rate coefficient depends on collisional (kinetic) energy and where the total
energy available to the products is deposited (i.e., into product translational energy or
product internal energy). The angular distribution of product molecules can also give
information about the nature of the reaction process. For example, if the A + B collision
forms a long-lived (i.e., on rotational time scales) collision complex, the product C
molecules display a very isotropic angular distribution.
In beam-bulb experiments, one is not able to gain as much detailed information
because one of the reactant molecules B is not constrained to be moving with a known
fixed velocity in a specified direction when the A + B → C collisions occur. Instead, the
B molecules collide with A molecules in a variety of orientations and with a distribution
of collision energies whose range depends on the Maxwell-Boltzmann distribution of
kinetic energies of the B molecules in the bulb. The advantage of beam-bulb experiments
is that one can achieve much higher collision densities than in crossed-beam experiments
because the density of B molecules inside the bulb is much higher than are the densities
achievable in a beam of B molecules.
There are cases in which the beam-bulb experiments can be used to determine how
the reaction rate depends on collision energy even though the molecules in the bulb have
51
a distribution of kinetic energies. That is, if the species in the beam have much higher
kinetic energies than most of the B molecules, then the A + B collision energy is
primarily determined by the beam energy. An example of this situation is provided by so-
called guided-ion beam experiments in which a beam of ions having well-specified
kinetic energy E impinges on molecules in a bulb having a temperature T for which kT
<<E. Fig. 8.6 illustrates data that can be extracted from such an experiment.
52
Figure 8.6 Collision-Induced Dissociation Data Showing Cross-Section as a Function of
Collision Energy.
In Fig. 8.6, we illustrate the cross-section σ (related to the bimolecular rate constant k by
σ v = k, where v is the relative collision speed) for production of Na+ ions when a beam
of Na+(uracil) complexes having energy E (the horizontal axis) collides with a bulb
containing Xe atoms at room temperature. In this case, the reaction is simply the
collision-induced dissociation (CID) process in which the complex undergoes
unimolecular decomposition after gaining internal energy in its collsions with Xe atoms:
Na+(uracil) → Na+ + uracil.
The primary knowledge gained in this CID experiment is the threshold energy E*; that is,
the minimum collision energy needed to effect dissociation of the Na+(uracil) complex.
This kind of data has proven to offer some of the most useful information about bond
dissociation energies of a wide variety of species. In addition, the magnitude of the
reaction cross-section σ as a function of collision energy is a valuable product of such
experiments. These kind of CID beam-bulb experiments offer one of the most powerful
and widely used means of determining such bond-rupture energies and reaction rate
constants.
C. Other Methods
53
Of course, not all chemical reactions occur so quickly that they require the use of
fast lasers to follow concentrations of reacting species or pump-probe techniques to
generate and probe these molecules. For slower chemical reactions, one can use other
methods for monitoring the relevant concentrations. These methods include
electrochemistry (where the redox potential is the specie's signature) and NMR
spectroscopy (where the chemical shifts of substituents are the signatures) both of whose
instrumental response times would be too slow for probing fast reactions.
In addition, when the reactions under study do not proceed to completion but exist in
equilibrium with a back reaction, alternative approaches can be used. The example
discussed in Chapter 5 is one such case. Let us briefly review it here and again consider
the reaction of an enzyme E and a substrate S to form the enzyme-substrate complex ES:
E + S ? ES.
In the perturbation-type experiments, the equilibrium concentrations of the species are
"shifted" by a small amount δ by application of the perturbation, so that
[ES] = [ES]eq -δ
[E] = [E]eq + δ
[S] = [S]eq + δ.
54
Subsequently, the following rate law will govern the time evolution of the concentration
change δ:
- dδ/dt = - kr ([ES]eq -δ) + kf ([E]eq + δ) ([S]eq + δ).
Assuming that δ is very small (so that the term involving δ2 cam be neglected) and using
the fact that the forward and reverse rates balance at equilibrium, this equation for the
time evolution of δ can be reduced to:
- dδ/dt = (kr + kf [S]eq + kf [Eeq]) δ.
So, the concentration deviations from equilibrium will return to equilibrium
exponentially with an effective rate coefficient that is equal to a sum of terms:
keff = kr + kf [S]eq + kf [Eeq].
So, by following the concentrations of the reactants or products as they return to
their equilibrium values, one can extract the effective rate coefficient keff. Doing this at a
variety of different initial equilibrium concentrations (e,g., [S]eq and [E]eq), and seeing
how keff changes, one can then determine both the forward and reverse rate constants.
