Subject: Chem Topic: Finding Transition Structures

From: David Young <[email-masked]>
Date: Wed, 31 Dec 1997 15:37:01 -0600 (CST)

Hello all,

	I have written the following short essay for my users.  I am 
posting it here because several of the subscribers to this list also
asked if I had any information on finding transition structures.
Feel free to send me any comments or suggestions.

	My compilation of chemical topics can be accessed via the web 
at URL http://www.auburn.edu/~youngd2/topics/contents.html

---------------------------------------------------------------------------

			Finding Transition Structures

   			          David Young

		        Division of University Computing
		  	         144 Parker Hall
		          Auburn University, AL  36849

INTRODUCTION

	A transition structure is the molecular species that is represented
by the top of the potential energy curve in a simple one dimensional 
reaction coordinate diagram.  The energy of this species is needed in order
to determine the energy barrier to reaction and thus the reaction rate.
The geometry of a transition structure is an important piece of information
for describing the reaction mechanism. 

	Short of determining an entire reaction coordinate, there are a 
number of structures and their energies that are important to defining a
reaction mechanism.  For the simplest single step reaction, there would be
five of these structures.

1.  The reactants separated by large distances.
2.  The van der Waals complex between the reactants.
3.  The transition structure.
4.  The van der Waals complex between the products.
5.  the products separated by large distances.

	A transition structure is mathematically defined as the geometry
which has a zero derivative of energy with respect to moving every one
of the nuclear coordinates and has a positive second derivative of energy
for all but one geometric movement which has a negative curvature.  
Unfortunately, this description describes many structures other than
a reaction transition such as an eclipsed conformation or the intermediate
point in a ring flip or any structure with a higher symmetry than the ground
state of the compound.

	Predicting what a transition structure will look like (without the
aid of a computer) is difficult for a number of reasons.  Such a prediction 
might be made based on a proposed mechanism which is incorrect.
The potential energy surface around the transition structure is often much
more flat than the surface around a stable geometry ... thus there may be
large differences in the transition structure geometry between two
seemingly very similar reactions and with very little differences in energy.

	Computationally it has been possible to determine transition 
structures for many years, although not always easy.  Experimentally
it has only recently become possible to examine reaction mechanisms directly 
using femtosecond pulsed laser spectroscopy.  It will be some time before 
these techniques can be applied to all of the compounds that are accessible
computationally.  Further more, these techniques yield vibrational information
rather than an actual geometry for the transition structure.

MOLECULAR MECHANICS PREDICTION

	Traditionally, molecular mechanics has not been the method of
choice for predicting transition structures.  However, since it is the only
method viable for many large molecules, some efforts have been made to
predict transition structures.  Since the bonds are explicitly defined
in molecular mechanics methods, it is not possible to simply find a point
that is an energy maximum.  The technique most often used (i.e. for a
atom transfer) is to first plot the energy curve due to stretching a bond
that is to be broken (without the new bond present), then plot the energy 
curve due to stretching a bond that is to be formed (without the old
bond present).  The transition structure is then defined as the point
at which these two curves cross.  Since molecular mechanics methods were
not designed to describe bond breaking and other reaction mechanisms, 
these methods are most reliable where a class of reactions has been tested
against experimental data to determine its applicability and perhaps
a suitable correction factor.

	The rest of the techniques mentioned in this document are applicable
to semiempirical, density functional theory (DFT) and ab initio techniques.

LEVEL OF THEORY

	As a general rule of thumb, transition structures are more difficult
to describe than equilibrium geometries.  As such, lower levels of theory
such as semiempirical methods, DFT using a local density approximation (LDA),
and ab initio methods with small basis sets do not generally describe
transition structures as accurately as they describe equilibrium geometries.
There are of course exceptions to this, but they must be identified on 
a case by case basis.  

	The best way to predict how well a given level of theory will 
describe a transition structure is to look up results for similar classes
of reactions.  Tables of such data are provided by Hehre in the reference
listed below.

USE OF SYMMETRY

	As mentioned above, a structure with a higher symmetry than is 
obtained for the ground state may satisfy the mathematical criteria 
defining a reaction structure.  In a few rare (but happy) cases, the
transition structure can be rigorously defined by the fact that it
should have a higher symmetry.  An example of this would be the symmetric
SN2 reaction

		F + CH3F -> FCH3 + F

In this case the transition structure must have D3h symmetry with the two
F atoms arranged axially and the H atoms equatorial.  In fact,
the transition structure is the lowest possible energy compound that 
satisfies this symmetry criteria.  

	In this case, the transition structure can be found by forcing
the structure to have the correct symmetry then optimizing the geometry.
This means geometry optimization rather than transition structure finding 
algorithms are used.  This is a benefit because geometry optimization 
algorithms are generally more stable and reliable than transition structure
algorithms.

	For systems where the transition structure is not defined by
symmetry it is often best to ensure that the starting geometry does not
have any symmetry.  This helps avoid converging to a solution which is
an energy maximum of some other type such as an eclipsed conformation.

OPTIMIZATION ALGORITHMS

	If a program is given a molecular structure and told to find 
a transition structure, it will first compute the Hessian matrix (the 
matrix of second derivatives of energy with respect to nuclear motion).
The nuclei are then moved in a manner which increases the energy in 
directions corresponding to negative values of the Hessian and decreases
energy where there are positive values of the Hessian.  This procedure has
several implications.

	This is a quasi-Newton technique which implicitly assumes that the 
potential energy surface has a quadratic shape.  Thus the optimization will 
only be able to find the correct geometry if the starting geometry is 
sufficiently close to the transition structure geometry to make this a 
valid assumption.  Quasi-Newton techniques are generally more sensitive
to the starting geometry than the synchronous transit methods discussed 
below.  One good way to get a structure close to the correct 
transition structure is to use a transition structure from a very similar 
system (i.e. the same reaction with different functional groups).

	Simplex optimizations have been tried in the past.  These do not
assume a quadratic surface, but require far more computer time and
are seldom incorporated in commercial software.  Due to the unavailability
of this method to most researchers, it will not be discussed further
here.

	The optimization of a transition structure will be much faster
using methods for which the Hessian can be analytically calculated.
For methods which incrementally compute the Hessian (i.e. the Berny algorithm)
it is advantageous to start with a Hessian from some simpler calculation, 
such as a semiempirical calculation.

FROM STARTING AND ENDING STRUCTURES

	Since transition structure calculations are so sensitive to the
starting geometry, a number of techniques for finding reasonable starting
geometries have been proposed.  One very useful technique is to start
from the reactant and product structures which are more easily obtained
than transition structures.

	The simplest way to guess the shape of a transition structure
is to assume that each atom is directly between the position where it 
starts and the position where it ends.  This linear motion approximation 
is called linear synchronous transit (LST).  This is a good first 
approximation, but it has its failings.  Consider the motion of an atom
which is changing bond angle with respect to the rest of the molecule.
The point half way between its starting and ending positions on the line
connecting those positions will give a shorter than expected bond length
and thus be (perhaps significantly) higher in energy.  

	The logical extension of this technique is the quadratic synchronous
transit method (QST).  These methods assume the coordinates of the atoms
in the transition structure will lie along a parabola connecting the 
reactant and product geometries.  QST generally gives some improvement
over LST although it may be a very slight improvement.

	Many programs allow the user to input a weighting factor (i.e. to 
give a structure that is 70% products and 30% reactants).  This allows the 
application of the Hammond postulate that the transition structure will look 
more like the reactants for an exothermic reaction and more like the products
for an endothermic reaction.

	These techniques have been very useful for simple reactions, but
have limitations.  The down side is that each of these, even at their
best is designed around the assumption that the reaction is a single
step with a concerted motion of all atoms.  For multi-step reactions,
these techniques can be used individually for each step.  For a reaction
which has only one transition structure but the motion is not concerted
(i.e. breaking one bond then forming another) it may be better to 
use starting geometries created by hand or eigenvalue-following.

	There are distinct differences in the way these methods are 
implemented in specific software packages.  Some software packages will
require the user to choose a transit method to obtain a starting geometry
then run a separate calculation with a quasi-Newton method.  Other software
packages will have an automated way of runing the transit method calculation
followed by a quasi-Newton calculation.  There have even been algorithms
proposed for allowing the program to make decisions concerning which
method to use at each step of the optimization. 

REACTION COORDINATE TECHNIQUES

	A transition structure is of course a maximum on the reaction
pathway.  One well defined reaction path is the least energy or intrinsic
reaction path (IRC).  Quasi-Newton methods oscillate around the IRC path
from one iteration to another and several groups have proposed methods
for obtaining the IRC path from the quasi-Newton optimization.  

	Likewise a transition structure can be obtained by following the 
reaction path from the equilibrium geometry to the transition structure.
This technique is known as eigenvalue-following because the user specifies
which vibrational mode should lead to a reaction given sufficient kinetic
energy.  This is not the best way to obtain an IRC, nor is it the fastest
or most reliable way to find a transition structure.  However, it has the 
advantage of not making assumptions about concerted motions of atoms
or what the transition structure will look like.

	Another technique is to use a pseudo reaction coordinate.  This
can be quite a bit of work for the user and requires more computer time
than most of the other techniques mentioned.  However it has the advantage
of being very reliable and thus will work when all other techniques have
failed.  A pseudo reaction coordinate is calculated by first choosing
a geometric parameter intimately involved in the reaction (such as the 
bond length for a bond that is being formed or broken).  A series of 
calculations is then run in which this parameter is held fixed at various
values from those in the reactants to those in the products and all
other geometric parameters are optimized.  This does not give a true
reaction coordinate but an approximation to it which matches the true
reaction coordinate perfectly only at the equilibrium geometries and 
transition structure.  Typically the highest energy calculation from this
set is used as the starting geometry for a quasi-Newton optimization.
In a few rare cases involving very flat potential surfaces the quasi-Newton
optimization may still fail.  In this case, the transition structure can
be calculated to any desired accuracy (within the theoretical model)
by finding the energy maximum by varying the chosen geometric parameter in
successively smaller increments.

POTENTIAL SURFACE SCAN

	The reaction coordinate is one specific path along the complete
potential energy surface associated with the nuclear positions.  It is
possible to do a series of calculations representing a grid of points on
the potential energy surface.  The saddle point can then be found
by inspection or more accurately by using mathematical techniques to
interpolate between the grid points.

	This type of calculation does reliably find a transition structure.
However, it requires far more computer time than any of the other techniques.
As such, this is really only done when the research requires obtaining
a potential energy surface for reasons other than just finding the
transition structure. 

SOLVENT EFFECTS

	It is well known that reaction rates can be affected by the
choice of solvent.  Solvent interactions can affect the energy of the
transition structure significantly and generally only slightly change
the transition structure geometry.  All of the techniques for finding 
transition structures can be used when solvent effects are being included
in the calculation.  The presence of solvent interactions does not change
the manner in which transition structures are found at all (although it
might change the results).

VERIFING THAT THE CORRECT GEOMETRY WAS OBTAINED

	The primary means of verifying a transition structure is to
compute the vibrational frequencies.  A saddle point should have one
negative frequency.  The vibrational motion with this negative frequency
is the motion going towards reactants in one direction and products
in the other direction.  

	It is also always important to look at the transition structure 
geometry to make sure that it is the reaction transition and not the 
transition in the middle of a ring flip or some other unintended process.
If it is not clear from the geometry, that the transition structure is 
correct, displaying an animation of the transition vibrational mode
should make it very clear.

	It is possible that a transition structure calculation will
give two negative frequencies (a second order saddle point) or more.
This gives a little bit of information about the potential energy surface
but it is extremely unlikely that such a structure has any significant
bearing on how the reaction occurs.  This type of structure will often
be found if the starting geometry was given a higher symmetry than the 
transition structure should have.

OBTAINING A REACTION RATE

	It is not the purpose of this document to give a detailed 
description of methods for obtaining reaction rates.  However, since
a reaction rate calculation is often the next step after finding a
transition structure, some of the issues involved will be mentioned here
for the sake of completeness.

	The simplest way to get a reaction rate is to use the activation
energy in the Arrehenius equation.  The preexponential factor can be
obtained from experimental observations or some simple theoretical 
method such as the kinetic theory of gasses.  To a first approximation
the activation energy can be obtained by subtracting the energies of the
reactants and transition structure.  A readily obtained additional 
correction to these energies is obtained by the addition of the zero point 
vibrational energy. 

	Simply using the activation energy assumes that the only way a 
reaction occurs is along the intrinsic reaction coordinate.  It would be
more correct to consider that reactions may occur which go through a 
geometry very similar to the transition structure as well.  Variational
transition state calculations take this into account.  These
calculations may require using the vibrational frequencies for the
transition structure, the entire reaction coordinate or the entire
potential energy surface.  These calculations can also take into account
tunneling through the reaction barrier.  These calculations can give
good results, but are very sensitive to subtle details like using a mass
weighted coordinate system to specify the geometry.

	Dynamical studies can be done to examine how the path and orientation
of approaching reactants affects the reaction rate.  These studies often
start with a potential energy surface which was obtained from ab initio
calculations.  The amount of work necessary to study a reaction with 
these techniques may be far more than the work done to get the potential
energy surface, which was not a trivial task in itself.

CHECKLIST OF METHODS FOR FINDING TRANSITION STRUCTURES

	Many techniques for finding transition structures are discussed
above.  The following is a listing of each of these starting with the
ones which are easiest to do and most often successful.  In other words,
start with number 1 and continue down the list until you find one that works.

1.  If the system can only feasibly be modeled by molecular mechanics
use the potential energy curve crossing technique.

2.  If the transition state can be defined by symmetry, do a normal 
geometry optimization calculation with the symmetry constrained.

3.  If you have the structure of the intermediate for a very similar
reaction, use that structure with a quasi-Newton optimization.

4.  Quadratic synchronous transit followed by quasi-Newton.

5.  Linear synchronous transit followed by quasi-Newton.

6.  Try quasi-Newton calculations starting from structures that look
like what you expect the transition structure to be like and have no
symmetry.  This is a skill which improves as you become more familiar with
the mechanisms involved, but requires some trial and error work even for 
the most experienced researchers.

7.  Eigenvalue-following.

8.  Pseudo reaction coordinate with one parameter constrained followed
by a quasi-Newton optimization.

9.  Pseudo reaction coordinate with one parameter constrained using
successively smaller steps for the constrained parameter until the 
desired accuracy is reached.

10.  Go back to options 8 & 9 and constrain a different parameter.

11.  Consider the fact that some reactions have no barrier.  You might also
be making incorrect assumptions about what the reaction mechanism should be.
Consider these possibilities and start over.

12.  Switch to a higher level of theory and start all over again.

13.  Obtain the transition structure from the entire potential energy surface.  
It is questionable if there can be any case where this is the only option
but it should work as a desperate last resort.

	Once you are experienced at finding transition structures for
a particular class of reactions, you will probably go directly to the 
technique that has been most reliable with those reactions.  Until that
time, this sequence is the authors best advice for finding a transition 
structure with the least amount of work for you and the computer.

REFERENCES

A good discussion of the issues involved and many tables of performance
data can be found in
W. J. Hehre "Practical Strategies for Electronic Structure Calculations"
Wavefunction (1995)

A nice discussion from the stand point of the potential energy surface
starts on page 240 of
A. R. Leach "Molecular Modelling Principles and Applications" Longman (1996)

For more information on synchronous transit methods see
C. Peng, H. B. Schlegel Israel Journal of Chemistry 33, 449 (1993)

Obtaining transition structures from molecular mechanics is discussed in 
F. Jensen J. Comp. Chem. 15, 1199 (1994)

---------------------------------------------------------------------------