Transition State Search (QST2 & QST3) and IRC with Gaussian09

Theoretical evaluation of a reaction mechanism is all about finding the right transition states (TS) but there are no guarantees within the available methods to actually find the one we need. Chemical intuition in the proposal of a mechanism is paramount. Let’s remember that a TS is a critical point on a Potential Energy Surface (PES) that is a minimum in every dimension but one. For a PES with more than two degrees of freedom, a hyper-surface, envisioning the location of a TS is a bit tricky, in the case of a three dimensional PES (two degrees of freedom) the saddle point constitutes the location of the TS as depicted in figure 1 by a section of a revolution hyperboloid.

400px-Saddle_point

Fig1. Saddle point on a surface (min in one direction; max in the other)

Fig 1a Pringles chips -Yuck-. They exhibit a maximum on the direction parallel to the screen and a minimum on the direction perpendicular to the screen at the same point.

Fig 1a Pringles chips -Yuck-. They exhibit a maximum on the direction parallel to the screen and a minimum on the direction perpendicular to the screen at the same point.

The following procedure considers gas phase calculations. Nevertheless, the use of the SCRF keyword activates the implicit solvent calculation of choice in order to evaluate to some degree the solvent influence on the reaction energetics at different temperatures with the use of the temperature keyword.

The first step consists of a high level optimization of all minimums involved, such as reagents, products and intermediates, with a subsequent frequency analysis that includes no imaginary eigenvalues.

In order to find the structures of the transition states we use in Gaussian the Synchronous Transit-guided Quasi-Newton method [1] through the keywords QST2 or QST3. In the former case, coordinates for the reagents and products are needed as input; for the latter keyword, coordinates for the TS structure guess is needed also.

QST2)

%chk=file.chk
%nprocshared=n
%mem=nGB

#p opt=(qst2,redundant) m062x/6-31++G(d,p) freq=noraman Temperature=373.15 SCRF=(Solvent=Water)

Title card for reagents

Q M
Cartesian Coordinates for reagents
blank line
Title card for products

Q M
Cartesian Coordinates for products
blank line

QST3)

%chk=file.chk
%nprocshared=n
%mem=nGB

#p opt=(qst3,redundant) m062x/6-31++G(d,p) freq=noraman Temperature=373.15 SCRF=(Solvent=Water)

Title Card for reagents

Q M
Cartesian Coordinates for reagents
blank line
Title card for products

Q M
Cartesian Coordinates for products
blank line–
Title card for TS
Q M
Cartesian Coordinates for TS
blank line

NOTE: It is fundamental that the numbering order is kept constant throughout the molecular specifications of all two, or three, input structures. Hence, I recommend to build a set of molecules, save their structure, and then modified the coordinates on the same file to produce the following structure, that way the number for every atom will remain the same for every step.

As I wrote above, there are no guarantees of finding the right TS so many attempts are probably needed. Once we have the optimized structures for all the species involved in our mechanistic proposal we can plot their energies very simply with MS Excel the way we’ve previously described in this blog (reblogged from eutactic.wordpress.com)

Once we’ve succeeded in finding the structure of our TS we may run an Internal Reaction Coordinate (IRC) calculation. This calculation will connect the TS structure to those of the products and the reagents. Initial constant forces are required and these are commonly retrieved from the TS calculation checkpoint file through the RCFC keyword.

%chk=QST3_2p.chk
%nprocshared=8

#p m062x/6-31++G(d,p) IRC=(Maxpoints=50,RCFC,phase=(2,1))Temperature=373.15 SCRF=(Solvent=Water) geom=allcheck

Title Card

Q M
blank line

Finally, the IRC path can be visualized with GaussView from the Results menu. A successful IRC will link both structures along a single reaction coordinate proving that both reagents and products are linked by the obtained TS.

Hat tip to Howard Diaz who has become quite skillful in calculating these mechanisms as proven by his recent poster at the XII RMFQT a couple of weeks back. And as usual thanks to everyone who reads, comments, likes, recommends, rates and shares my silly little posts.

About joaquinbarroso

Theoretical chemist in his early thirties, in love with life and deeply in love with his woman. I love science, baseball, literature, movies (perhaps even in that order). I'm passionate about food and lately wines have become a major hobby. In a nutshell I'm filled with regrets but also with hope, and that is called "living".

Posted on November 27, 2013, in Computational Chemistry, Gaussian, Models, Research, Theoretical Chemistry, White papers and tagged , , , , , , , , , , , . Bookmark the permalink. 20 Comments.

  1. Anya Skitnevskaya

    It’s a little misprint. In QST3 input you put: opt=(qst2…
    Nevertheless. Already for a long time I wanted to say thank you for your blog, you do a great and important work! :)

  2. I am trying to get a transition state having a complex molecular system with 4 H-bond. the The groups invoved for interaction are C=O and NH2. The distance between them is 4.5 angstrom initially. I got a transition state from it with 48 kCalmol(-). It seems that I have missed a transition state in between. What procedure should I take to get the transition state in between the H bonded initial complex and the proton shuttling TS?

  3. Thanks in support of sharing such a pleasant idea, paragraph is
    good, thats why i have read it fully

  4. Hola, estoy intentado proponer un estado de transición que es básicamente una disociación para el (CH3)3CC(O)OO—-NO2 pero cuando obtengo el TS la frecuencia negativa que encuentro es muy pequeña así que necesito mejorar esa frecuencia. Los keywords que estoy usando son
    # opt=(calcfc,tight,ts,noeigentest) freq ub3lyp/6-31++g(d,p)
    Me podría dar algunas recomendaciones que me ayuden a mejorar el cálculo?
    Muchas gracias.

    • Hola Diana!

      A que te refieres con muy pequeña? Ni en valor ni en intensidad (IR o RAMAN) necesitas que sea “grande” si corresponde al estado vibaracional que conecte tus reactivos con tus productos.
      Por lo que veo tu nivel de teoría es adecuado, si bien el uso de b3lyp pudiera acaso ser cuestionable aunque nosotros hemos obtenidos buenos resultados (reproducciones experimentales) con dicho funcional.

      Saludos y estoy al pendiente

      • La frecuencia IR es -54.23, lo que quiero observar es la separación entre el O—N del sistema así que digamos que no necesariamente que sea grande pero que sea la frecuencia para este vibración. Y bueno precisamente estoy buscando este estado de transición porque quiero compararlo con lo obtenido experimentalmente. No se si queda claro?

      • Pues si esa frecuencia de -54.23 corresponde a la separación O-N que buscas pues ya terminaste! O estoy entendiendo mal?

  5. Excellent article. Keep writing such kind of information on your blog.
    Im really impressed by it.
    Hello there, You’ve performed an incredible job. I will
    definitely digg it and in my opinion recommend to my friends.
    I am confident they will be benefited from this website.

  6. Buenas Joaquin:
    Encontre un TS a nivel HF/3-21G* con una frecuencia imaginaria de -186 cm-1 y el calculo del IRC me da perfecto. Sin embargo al querer reoptimizar esa estructura usando M062X/6-311+G** me genera una estructura con una frecuencia negativa muy pequeña -39 cm-1 y el IRC solo me da dos pasos y se corta. No estoy segura si en realidad es un TS.

  7. ‘m working in IRC calculation in Gaussian 09. The job was terminated with the following error message. How can I resolve it?

    The error message: Maximum number of corrector steps exceded.
    Error termination via Lnk1e in g09/l123.exe

  8. Hi
    If it possible one exaple about QST2 and QST3 and RIC

    Thanks a lot

  9. Hi thanks for the useful information.

    I am trying to model a complex reaction and found the optimised reagents and products (from opt jobs) to be higher in energy than the end structures from IRC curves. I think this is a rather odd problem as I think they should be lower in energy if anything. Is there an immediate reason why? Freq jobs show no imaginary frequencies for the optimised reagents/products, and indeed, the opt job went to completion – but I understand they might have been optimised at a higher minimum.

    The transition states were aquired with opt=qst3 and show only one imaginary frequency. The IRCs do show reagent/product formation for reverse and forward directions, but of course the bond lengths and angles aren’t exact which is a little expected. I’ve been using the B3LYP/6-311++G(d,p) level as this showed the least sharing of MOs due to lack of diffuse functions on lesser basis sets.

    The root line for the optimised reagents is:
    #n B3LYP/6-311++G(d,p) Opt Freq

    The root line for the transition states is:
    #n B3LYP/6-311++G(d,p) Opt=QST3 Freq

    The root line for the IRCs is:
    # B3LYP/6-311++G(d,p) irc=(CalcFC,forward,maxpoints=79,stepsize=10) geom=check

    I used 79 max points as this was the maximum path steps found for the run when I initially set 200 maxpoints.

    Any help would be most appreciated,
    Leighton

    • Hi Leighton

      We have come across this problem too. It is probably due to a complex potential energy surface. Unfortunately the only advice I can give you is to play -systematically- a bit with geometries and levels of theory. Using the NoEigenTest keyword might help you but it will consume more computing time.

      I hope this helps, let me know!

      • Hi Joaquin,

        Unfortunately, the NoEigenTest did not change the situtation – although the good news is that the NoEigenTest ran fine and jobs went to completion! Thanks for the suggestion all the more.

        So now another question along similar lines if I may…I would like to calculate the energy surface, starting with the saddle point and then moving down the IRCs steepest descent forwards and reverse to the local minima (probably done over many jobs). How do I do this? I can only find Scan keywords where the only examples provided are dihedral rotations, but not the generation of a saddle point surface. Can you recommend a keyword in Gaussian 09 that can achieve this please?

        I would like to specifically refer to the type of surface by referencing Figure 1 in the following freely accessible article: http://chem.wayne.edu/schlegel/Pub_folder/265.pdf

        Best regards,
        Leighton

  10. Kudos

    I use your page for tutoring students…

    Manythanks…

    k

  11. Hi Dr. Joaquin,
    Thanks a lot for your useful blog.
    please, I want to know how can I convert temperature and pressure in G03W, we know that the default temperature is 298.15K and pressure 1 atm. I need to make my calculations at another conditions (e. g. T=0K). my calculations are not frequencies or thermochemistry properties, they are energies for a system.

    I used
    # … Temperature=0K

    but it is failed

    • Hello,

      The idea of pressure and temperature for a single molecule in a vacuum lacks any physical meaning. The Temperature keyword is used for the thermochemistry step but it cannot be used for the energy of a single molecule.

      I hope this helps

  12. Hello,

    I’ve been trying to find the TS barrier energies associated with a linear N—-O bond breakages. I’ve been using TS with noeigentest and have had one of two results:

    A TS showing the suspected product molecule (bond unbroken), typically with neg freqs (1-3) showing rotational degrees of freedom.

    A TS showing an infinitely separated complex, usually ending in many errors.

    TS search is a frustrating (and rewarding) venture to say the least. Nucleophilic substitutions are pretty simple to find, but bond breakage has been giving me troubles. Any tips for moving forward? I haven’t had much luck with typical hybrid density functionals (b3lyp) with reasonably sized basis sets.

    Thanks for any advice,

    Daniel

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 1,554 other followers

%d bloggers like this: