204

Straub

where ε‡PR is the activation energy required to pass from the reactant state to the product state. The angular frequency of the reactant well is ω0.

C. Computing the TST Rate Constant

What knowledge of the system is necessary to compute kRPTST?

1.We need to have a good estimate of the energy of the system as a function of the positions and momenta of all atoms in the system [13,14].

2.It is necessary to compute the average over the phase space of the system.

3.We must be able to define the reaction coordinate along which the transition state theory dividing surface is defined.

Each of these requirements can be difficult to meet for a complex biomolecular system. Each of these points is addressed in this chapter.

A variety of methods for finding reaction paths in simple chemical systems have been proposed. Good review articles summarizing those methods can be found [8,15,16]. An excellent historical overview of these methods is provided by Anderson [17]. Here we focus our discussion on those methods that have had the widest application to largescale biomolecular systems and that hold the greatest promise for further development.

III. CORRECTIONS TO TRANSITION STATE THEORY

The assumptions of transition state theory allow for the derivation of a kinetic rate constant from equilibrium properties of the system. That seems almost too good to be true. In fact, it sometimes is [8,18–21]. Violations of the assumptions of TST do occur. In those cases, a more detailed description of the system dynamics is necessary for the accurate estimate of the kinetic rate constant. Keck [22] first demonstrated how molecular dynamics could be combined with transition state theory to evaluate the reaction rate constant (see also Ref. 17). In this section, an attempt is made to explain the essence of these dynamic corrections to TST.

Transition state theory assumes that once the system reaches the transition state the system dynamics will carry the activated reactant to product, where it will be deactivated. It assumes that the process of converting activated reactants to products is perfectly efficient. Thinking about the system dynamics, we understand that that will not always be the case. For a reaction system in an environment in which the collision rate (or friction) is low, the activated system may cross into the product state and recross the transition state surface back to the reactant state many times before undergoing collisions, losing energy, and becoming deactivated. Alternatively, when the collision rate (or friction) is very high, the activated system may be kicked back and forth across the transition state surface many times before being deactivated. Dynamics typical of both regimes are depicted in Figure 3. These dynamic processes in the low and high friction regimes can be effectively studied by using molecular dynamics simulations.

Either of the mechanisms of recrossing leads to inefficiency in converting reactant to product. How does this affect the reaction rate constant? Fewer activated reactants form

Reaction Rates and Transition Pathways

205

Figure 3 Dynamic recrossings in the low and high friction regimes. Recrossings back to the reactive state lead to a lowering of the rate constant below the transition state theory value.

products, so the rate constant will be lower than the TST estimate. This is summarized in the formula for the actual rate constant,

kexact κkTST

(12)

where κ is the transmission coefficient—a positive number less than or equal to 1.

A. Computing Using the Reactive Flux Method

In practice, we can compute κ as follows [19,23]. We start with a set of trajectories at the transition state q q‡. The momenta have initial conditions distributed according to the normalized distribution functions

206	Straub
P( )(Γ, γ) δ(q q‡)θ( p)p exp( β )	(13)

where Γ and γ are the nonreactive and reactive phase space degrees of freedom, respectively. In one set, P( )(Γ, γ), the trajectories initially have positive momenta (and at first move into the product well). In the complementary set, P( ) (Γ, γ), the trajectories initially have negative momenta (and at first move into the reactant well).

Using these distribution functions, we can write the reactive flux correlation function in the compact form

ˆ		( )	( )	‡	]	(14)
k(t) ∫dΓ dγ[P			(Γ, γ) P	(Γ, γ)]θ[q(t) q
ˆ						‡	]
Note that k(t 0)		1. What does this function measure? The function θ[q(t) q

follows each trajectory and counts 1 if the trajectory is in the product well and 0 otherwise. The calculation begins at t 0 with a number of trajectories distributed according to P( )(Γ, γ) and an equal number according to P( )(Γ, γ). The trajectories are followed

in time until they are deactivated in the reactant or product well. This is illustrated in Figure 4. Initially there may be rapid recrossings of the transition state, and this can lead

ˆ

to a rapid decay or ringing in k(t). After a time, all of the transient decay will have passed and only motion on the longest time scale—the time scale for activated barrier crossing—

ˆ

will be active. Eventually, k(t) will decay to zero on that long time scale. However, at the intermediate time scale, longer than the time scale for transient decay and shorter than

ˆ

the time scale for activated barrier crossing, the function k(t) will equal the transmission coefficient κ. Once κ is known, the total rate constant can be computed by multiplying κ by the TST rate constant.

If the assumptions underlying the TST are satisfied, the trajectories with initially

positive momenta will be trapped in the product well and those with initially negative

ˆ momentum will be trapped in the reactant well. That will result in a value of k(t) 1

and the rate constant k kTST.

If there are recrossings of the transition state, this will cause the positive contribution

ˆ

to k(t) to be somewhat less than 1 and the negative contribution to be somewhat greater

Figure 4 Reactive flux calculation for dynamics at low friction.

Reaction Rates and Transition Pathways		207
ˆ	TST	less than the transition
than 0, leading to a decay in k(t) 1 and a rate constant k κk
state theory estimate.

B. How Dynamic Recrossings Lower the Rate Constant

Consider a symmetrical double well (but the argument is easily generalized to the asymmetrical case). The evolution of a set of trajectories is schematically shown in Figure 5. In TST all trajectories crossing the transition state surface from reactants to products contribute to the reactive flux. However, when there are recrossings this is not the case. Some trajectories that are part of the ensemble of forward-moving trajectories actually originated as products. The contribution of those trajectories must be subtracted. Other trajectories started as reactant but will be deactivated as reactant. Those trajectories should not count at all. Here is how we can perform the counting that is done when computing the reactive

ˆ

flux k(t).

We assume that when the activated reactants cross the transition state a fraction P are deactivated as product and the remaining fraction 1 P recross the transition state surface [8,24]. If each fraction has roughly the same distribution of momenta as the original fraction, we can say that of the fraction 1 P that recross, P(1 P) will be deactivated in the reactant well and the remaining (1 P)2 will recross the transition state into the

Figure 5 The transition state ensemble is the set of trajectories that are crossing the transition state from reactant to product at equilibrium (shown as black dots). There are four types of trajectories, shown top to bottom in the diagram. (1) Starting as reactant, the trajectory crosses and recrosses the transition state and is deactivated as reactant. It does not add to the reactive flux. (2) Starting as reactant, the trajectory is deactivated as product. It adds 1 to the reactive flux. (3) Starting as product, the trajectory crosses and recrosses the transition state and is deactivated as product. Such a trajectory must be subtracted from the ensemble, so it counts 1 to the reactive flux. (4) Starting as product, the trajectory is deactivated as reactant. It does not contribute to the reactive flux.

208

Straub

product well, where a fraction P(1 P)2 will be deactivated as product. And so on. Adding up all contributions we find that a total fraction

	P P(1 P)2 P(1 P)4	1		(15)
		2 P

is deactivated as product.

But this is not the whole story! We not only need to know that a trajectory that crosses the transition state surface is eventually deactivated as product, we also need to know whether it originated from the reactant well! A trajectory that originates from the product well and ends up as product won’t contribute to the forward rate of reaction. Some of the trajectories did originate as product. We need to find that fraction and subtract it.

What is that negative contribution? We can follow the trajectories backward in time to find the well from which they originated. Of the number of trajectories initially moving from product to reactant, a fraction P is deactivated as reactant and a fraction 1 P recross the TST due to inertial motion or frequent collisions. A fraction P(1 P) will then be deactivated as product, and the remaining (1 P)2 will recross. And so on. The total fraction that is deactivated as product is

P(1 P) P(1 P)3	1	P	(16)
		P
2

Now we can compute the transmission coefficient [17,24]. It will be the difference between the positive and negative contributions, or

κ	1			1	P		P		(17)
	2 P						2 P
			2		P

Note that when P 1 we find that the assumptions of TST are met and κ 1. As the number of recrossings of the transition state increases, both P and κ decrease.

C. An Efficient Method for Computing Small Values of

In the very high and very low friction regimes, it might be that most trajectories do recross the transition state. It may also be that it takes a very long time to follow the system dynamics to the plateau region where κ can be measured. In such cases, it is also possible to compute the transmission coefficient in an approximate but accurate manner [24]. By placing an ‘‘absorbing boundary’’ at the transition state and simply following the trajectories until they recross or are deactivated, it is possible to estimate the value of κ with much less computational effort. Trajectories are started in the normal way as in the calculation of the reactive flux. However, when a trajectory recrosses the transition state, the run is stopped. By computing the fraction of trajectories (P ) that recross the transition state from the P( )(Γ, γ) distribution and the fraction (P ) that recross from the P( )(Γ, γ) distribution, we can estimate [24]

κ

P P

(18)

P P P P

When there are many recrossings, P and P may be much less than 1, few trajectories are integrated for long times, and the computational saving can be great. When the poten-