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.


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.



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

Title card for reagents

Cartesian Coordinates for reagents
blank line
Title card for products

Cartesian Coordinates for products
blank line



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

Title Card for reagents

Cartesian Coordinates for reagents
blank line
Title card for products

Cartesian Coordinates for products
blank line–
Title card for TS
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

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.


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

Title Card

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 late thirties, in love with life and deeply in love with his woman and children. 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. 31 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,

    • 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:

        Best regards,

  10. Kudos

    I use your page for tutoring students…



  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,


    • Hello Daniel,

      “S search is a frustrating (and rewarding) venture to say the least” you said it! It is indeed a bit of an art form.
      It seems you are on the right track but maybe changing functionals would go a long way. Try using M06-2x instead of b3lyp, that way you’ll get the long range interaction right. Also try the QSTn methods at least to a first approximation to the geometry of your TS.

      I hope this helps. Have a nice day!

  13. Hi, I’ve tried to run this calculation and i got this error ” Maximum number of corrector steps exceded.”.
    Could you help me?

  14. Hi DrJoaquin
    what can be done if an IRC run in g09 goes in the same direction for both the reverse and forward keywords?I try with keyword phase but see the
    same problem.what time does use keyword phase? which had atoms chosen? I try with atoms that formed the new bond between two molecules. is it correct? please guide me.

  15. Hi Dr Joaquin many thanks for your good answers,please if you illustrate qst2 and irc with clear ep.

  16. Hola, Joaquín. En primer lugar felicitarte por el blog, pero tengo algún problema con el cálculo. Estoy usando la siguiente línia de comandos:

    #P b3lyp/6-31+G** Opt=(qst2,redundant) geom=allcheck gfprint gfinput gfoldprint SCRF=(PCM,SOLVENT=WATER)

    Pero me da el siguiente error:

    Number 0
    Base 40448
    End 65536
    End1 65536
    Wr Pntr 40448
    Rd Pntr 40448
    Length 25088
    Error termination in NtrErr:
    NtrErr Called from FileIO.

    He leído por ahí que puede ser que falte memoria de cálculo, pero estoy usando 33 GB. Además, no me funciona con otras bases, ni tan siquiera con una hf/3-21G.
    Tienes alguna idea de dónde puede estar el problema?

    Muchas gracias

    • Posiblemente se trata de un error en la geometría, en tu línea de entrada dice que estás leyendo la geometría del archivo chk, eso es correcto?
      La última línea del error indica que estás leyendo desde un archivo que no necesariamente existe.
      Si en tu carro estás definiendo la creo metría de entrada elimina el comando geom=allcheck y con eso debería funcionar

  17. Hola joaquin, buenas noches, felicitaciones por el blog, me ha sido util.
    Estoy realizando un calculo IRC, cuando el proceso termina, me aparece este mensaje: this type of calculation cannot be archived. A que se debe esto?

  18. Dear Dr. Joaquin
    How are you doing? 🙂 It’s me, again. I have a question regarding the study of TS using the input “scan”. I was wondering if there’s a difference between doing a scan in gaseous phase and then optimizing the TS in a solvent and doing the scan in the solvent from the beginning.
    I run a calculation using both ways and got a small difference in the data, but I’m not sure it’s a significant difference. I hope you can illustrate me in the matter 🙂
    Thank you ever so much,
    Maria Fernanda (Caracas – UCV).

    • Hello María Fernanda,

      The problem I could find with merging both schemes is the cavity. When you optimize anything with a PCM cavity then you constrain the electron density (a main portion of it, anyway) to always be inside that cavity; If during the TS optimization some of it falls out then the program would stop due to an erroneous calculation of the solvation part, not necessarily that of the TS search.
      If you were able to get to significantly similar results in both approaches I would advice you to cosider that is not always the case.
      Of course the more complete method would be to use the solvent from the beginning but if for some reason PCM part crashes then I would settle for the second approach.

      I hope this helps. Have a nice day!

  19. Dear Dr. Joaquin,
    Thank-you so much for your blog, its helped me a lot throughout my Ph.D. and the supervision of other students.

    Recently I’ve been attempting to do an IRC calculation with no success. My system consists of a protein active site and as such, i have used modredundant to fix some coordinates. However, I appear to be encountering some errors in my calculations.

    This is the error that i get:
    ” Charge and Multiplicity card seems defective:
    Wanted an integer as input.”

    Which doesn’t make sense to me as my charge and multiplicity are inputted correctly, therefore, I’m wondering if this has something to do with using modredundant in an IRC calculation.

    What are your thoughts on this?

    Many thanks,

Leave a Reply

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

You are commenting using your 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

%d bloggers like this: