Subject: Summary: Loacting a transition state

From: [email-masked]
Date: Wed, 15 May 2002 16:21:47 -0500

Hi all,

I would like to thank those who responded my questions on locating a 
transition state.  All of them are quite instructive for me. Special 
thanks go to  David Gallagher and James J.P. Stewart who actually located 
the transition state for my system (See below). I could also locate a 
transition state.

Responses exhibited wide spectra of perception of the transition state 
search:

|---- Easy                                              David Gallagher, 
James J.P. Stewart, David Close
|
|---- Accurate methods                    Valentine Ananikov, Alexander 
Martins
|
|---- Difficult                                          Irene Newhouse, 
Bin Shan
|
|---- Black art                                      Markus Dittrich
|
|----More than one TS                     Stephen Bowlus
|
|---- Existence of TS                        David Shobe

 David Gallagher did find the transition state for the hydrolysis reaction 
of acetylcholine with water using PM5 method in CAChe version of MOPAC 
2002. Using PM3 in MOPAC, James Stewart proved that the structure derived 
> from the above PM5 method has only one negative frequency. I used the PM5 
TS structure with PM3 in Gaussin98 [#P RHF/PM3 opt=(TS,CalcFC)], however, 
the job stopped yielding two negative frequencies (Optimization stopped. 
-- Wrong number of Negative eigenvalues: Desired=  1 Actual=  2 ). The 
suggestion by Valentine Ananikov, Alexander Martins and Bin Shan 
[opt=(modredundant,TS,noeigen)] solved this problem. 

In the above, we assume that there is a transition state for a given 
reaction. David Shobe raised a question whether or not the transition 
structure exists, pointing out a no barrier reaction for TiCl4 + H2O --> TiCl4(H2O). The existence of a TS may be considered from the free energies for the 
reactants and the products. If the products are higher in free energy than 
the reactants, there must be a barrier. However, for the opposite case, 
there is no grantee that the barrier exists. Are anyone aware of any 
criteria for the existence of the transition state?

Regards,

Genzo

My posting was

Hi, 

I am trying to locate a transition state for a system about 30 atoms using 
an empirical
method. So far I haven't had any success. I would appreciate if you could
give me any suggestions. 

I am interested in a hydrolysis reaction pathway, using PM3 in Gaussian 
98: 

R-C-O-C-R' + H*-O'H ---> R-C-O-H*  + H-O'-C-R' 

I have done the following: 
1) opt=Ts ; This always failed due to incorrect number of negative
frequencies. 
2) Scan  ; I have not tried much. 
3) opt=ModRedundant: Scanning succeeded, but the population analysis after
the scanning seemed failed. The output did not indicate why it failed, but
the heading of population analysis appeared, followed by signal 11 error. 
4) opt=qst2; There is (are) one or two statement(s) "stationary point was
found" (the distances were the same for all, but the angles and torsions 
were slightly different in the two stationary points) in the output, but 
the optimization (?) kept going and failed; Optimization aborted
--- No acceptable step or Inconsistency: ModMin=   2. 

I optimized two reactant molecules individually and moved them close to 
make
a reaction. The products were also optimized. The order of the atoms are
exactly the same. 

Thanks, 

Genzo Tanaka 

Here are the responses for my questions

***
This should be easy, but what exactly is R-C-O-C-R meant to be?

An ether (R-CH2-0-CH2-R) or an anhydride (R-CO-O-CO-R) or what?

If you send me the exact structures, I'll try and run it for you.

David Gallagher   <[email-masked]>
***

Although, the old MOPAC PM3 hamiltonian seems to have difficulty in 
locating a transition state for your reaction, the latest PM5 method in 
MOPAC 2002 nailed it first go with a saddle calculation.

The attached Zip file includes pictures and the MOPAC output of the gas 
phase transition state for your hydrolysis reaction. In the picture, the 
bond thickness indicates the relative bond order and the numbers are the 
atom distances in Angstroms. The force calculation shows a single negative 
vibration and the IRCs verify that this is the only transition state along 
this particular reaction path for the concerted reaction. The new PM5 
method has been re-parameterized for higher accuracy and MOPAC 2002 
includes new algorithms. It seems to be much faster and more reliable for 
finding transition states than the older MOPAC methods such as PM3 and 
AM1. Also, I was able to model the transition state with the COSMO water 
solvent field. As expected, the energy was much lower.

Although, the old PM3 seems to work fine for a simple ester hydrolysis, 
the charged quaternary ammonium group seems to confuse it. The best 
transition state structure I got out of PM3 had a second negative 
vibration due to a methyl rotation. However, the bond orders, atom 
distances and negative vibration all indicated a proton transfer only (not 
the C-O bond break), so I don't think that this was the correct transition 
state for the reaction. 

If you need more information about MOPAC 2002 you can view the manual at http://www.cachesoftware.com/techsupport/mopac/ The pictures and calculations were done with the CAChe version of MOPAC 
2002.

I hope this helps
Regards
David Gallagher, Fujitsu   <[email-masked]>
***
Attached are two MOPAC data sets.  In the first, bug.mop, I used David's 
geometry from his PM5 result.  I changed PM5 to PM3, and ran the TS.  This 

runs without problem.

 From the ARC file, I made bug1.mop.  This runs the FORCE calculation and 
confirms that there is exactly 1 imaginary frequency.

So, in summary, both PM3 and PM5 (and, I assume, all the other methods) 
give more-or-less the same transition state.

As far as the system being difficult, my apologies, but what was 
difficult?

(I use MOPAC as stand alone, because I'm more comfortable with that.)

Best wishes,

Jimmy   <[email-masked]>
***

I'm not really an expert on TS optimizations, although I was willing to 
help
someone who was even further back on the learning curve than I was.  And
I've never had much luck with opt=qst2.  You've checked that the atoms
correspond properly between reactant and product, which is good.  I think
sometimes opt=qst2 has trouble if the path from reactants to products is
indirect: opt=qst2 will try to find a path that takes the reagent 
*through*
the other molecule, instead of going *around* it. 
 
I don't even know what signal 11 error is.  It sounds like a Unix message.
Two possibilities that come to mind are that there was a time limit for 
the
job which ran out, or that the system administrator "killed" it for some
reason.
 
Given what you've told me, I would try the opt=modredundant again, after I
found out what Signal 11 error means and had done whatever was necessary 
to
avoid it.  Oh, and you may want to send your question to the CCL list, 
since
someone else may have a better idea!
 
--Dave S.   <[email-masked]>
***
I think of these as being multi-step reactions.  Maybe you're trying to 
find
a TS that doesn't exist, because you skipped a step? 

> I am interested in a hydrolysis reaction pathway, using PM3 in 
> Gaussian 98:
>
> R-C-O-C-R' + H*-O'H ---> R-C-O-H* + H-O'-C-R'
>

It turns out that in one of the examples that poisoned my mind against
opt=qst2, the transition state was not found for an equally good reason: 
the
reaction was an association reaction TiCl4 + H2O --> TiCl4(H2O), and as it
turns out the energy is monotonic from reactants to product: there is no
activation barrier.

Anyway, I wish you luck in solving your transition state puzzles.

--David Shobe   <[email-masked]>
***
Transition states are always a difficult problem.  You will have to be 
pretty close to the actual state to get it to optimize.  Yes, that's as 
circular as it sounds!

I would strongly urge you to search the literature for similar transition 
states & use those structures to guide you in your initial guess.

Irene Newhouse   <[email-masked]>

***
The results of the first step are very common.  But this is not a problem.
You must look at the output with software that will show the vibration
frequencies.  Only one of them will involve motion along the coordinates
of the reaction path you are trying to follow.  Look at the other negative
frequencies.  You will see that they have nothing to do with your reaction
coordinates.  It may be some low frequency involving the wagging of an
exo-cyclic group.  Find the atoms involved in the unwanted motion.  Go to
the Standard Orientation near the end of the output and disturb these 
coordinates a little (maybe change them 0.1 A).  Run the OPT (TS) again
with these new coordinates.  You will quickly come down off the local
minimia you are stuck on and be at the true TS.
  Regards, Dave Close.   <[email-masked]>
***

First, I would further research whether the (apparently) four center 
reaction that is being done here is really feasible (or is that the 
question you are trying to answer?)  I would guess there are two TS that 
are encountered in this reaction.  For example, addition of water to form 
RCO-  ... +H2O'CR' (tight) ion pair, followed by hydrogen transfer to form 
the two neutral species.
 
Regardless of the mechanism, I would use something like SADDLE in Mopac or 
CHAIN in Ampac to get an estimate of the possible TS, then refine the TS 
using eigenvector following or some gentle gradient reduction method.  The 
general idea is to come up with an approximate starting geometry that is 
somewhere on the reaction coordinate (not just two molecules in space; 
this is where some inference about the probable mechanism is needed); a 
final product geometry that is on the other side of the hump, and then 
work both ends against the middle to find an approximate TS.  Then refine. 
 Then do IRC calculations to assure that you are actually doing the 
chemistry you think.  There is a section in the Mopac manual about all 
this, and it is included in an early JCAMD (around 1996, I think).
 
Sorry if this is redundant to what you know ... 
 
Stephen Bowlus    <[email-masked]>
***
I can't tell you much concrete. Finding transition states aka
looking for saddle point is a black art. I am currently
doing something similar for a 20 atom system using ab initio methods and
don't have much experience with semi-empirical ones.

In general you have to start as close as possible to the true
transitions state when starting a saddle point search. Otherwise
you won't be able to get it. You could make an educated guess or
first constrain some degrees of freedom you think are important
while looking for saddle points. Comparing the energies can give
you a hint in which direction you have to go.
Just moving two optimized molecules close to each other won't
do in 99.9% of the cases. Your transition state will be
'something inbetween'. Maybe a partial transfer of your groups,
partial bond formation... You will have to play around.
If you have a good guess about a possible reaction coordinate
contrain you system such that you mimick the reaction, monitoring
the energies. That should allow you to position your system as
close to the transition state as possible.

In any case. Once you think you found the transition state you
have to make sure that you have found a saddle point
by recalculating your hessian matrix which has to have one and
only one imaginary frequency. Otherwise you're back to where
you were before.

Good luck,

Markus Dittrich   <[email-masked]>
***
opt(TS, noeigentest) should help.

I suggest you HF/3-21G instead of PM3 for preliminary calculations
and B3LYP/6-31G* for more accurate data.

regards, Valentin Ananikov.   <[email-masked]>
***
Hi,

        First of all, it's a good idea use DFT methods instead of 
semiempirical ones. The results ara more reliable.

        You can use the keyword OPT=(TS, NoEigenTest) with a good guess 
for the geometry. Futhermore, the keyword OPT=QST3 it's a good option 
too. Don't forget use a good basis set, like 6-31G.

        Try these sugestions.

        Good luck,

        Alexander Martins.   <[email-masked]>
***
Hi,
    Genzo,try something like this

# opt(modredundant,TS,noeigen,CalcFC) rb3lyp/6-31g(d) nosymm

CalcFC calculates the initial force constants be computed at the first 
point 
and helps to determine which direction the atoms are to to move.The 
keyword 
noeigen and nosymm are also important.
   And it is also extremely important to specify a good initial 
guess.Finding transition state is by no means an easy task...according to 
my 
experience.
   Good luck.

Bin Shan   <[email-masked]>
***