The ‘art’ of finding Transition States Part 1

Guillermo CaballeroGuillermo Caballero, a graduate student from this lab, has written this two-part post on the nuances to be considered when searching for transition states in the theoretical assessment of reaction mechanisms. He’s been quite successful in getting beautiful energy profiles for organic reaction mechanisms, some of which have even explained why some reactions do not occur! A paper in Tetrahedron has just been accepted but we’ll talk about it in another post. I wanted Guillermo to share his insight into this hard practice of computational chemistry so he wrote the following post. Enjoy!

 

Yes, finding a transition state (TS) can be one of the most challenging tasks in computational chemistry, it requires both a good choice of keywords in your route section and all of your chemical intuition as well. Herein I give you some good tricks when you have to find a transition state using Gaussian 09 Rev. D1

METHOD 1. The first option you should try is to use the opt=qst2 keyword. With this method you provide the structures of your reagents and your products, then the program uses the quadratic synchronous transit algorithm to find a possible transition state structure and then optimize it to a first order saddle point. Here is an example of the input file.

link 0
--blank line--
#p b3lyp/6-31G(d,p) opt=qst2 geom=connectivity freq=noraman
--blank line--
Charge Multiplicity
Coordinates of reagents
--blank line--
Charge Multiplicity
Coordinates of products
--blank line---

It is mandatory that the numbering must be the same in the reagents and the products otherwise the calculation will crash. To verify that the label for a given atom is the same in reagents and products you can go to Edit, then Connection. This opens a new window were you can manually modify the numbering scheme. I suggest you to work in a split window in gaussview so you can see at the same time your reagents and products.

The keyword freq=noraman is used to calculate the frequencies for your optimized structure, it is important because for a TS you must only observe one imaginary frequency, if not, then that is not a TS and you have to use another method. It also occurs that despite you find a first order saddle point, the imaginary frequency does not correspond to the bond forming or bond breaking in your TS, thus, you should use another method. I will give you advice later in the text for when this happens. When you use the noraman in this keyword you are not calculating the Raman frequencies, which for the purpose of a TS is unnecessary and saves computing time. Frequency analysis MUST be performed AT THE VERY SAME LEVEL OF THEORY at which the optimization is performed.

The main advantage for using the qst2 option is that if your calculation is going to crash, it generally crashes at the beginning, in the moment of guessing your transition state structure. Once the program have a guess, it starts the optimization. I suggest you to ask the algorithm to calculate the force constants once, this generally improves on the convergence, it will take slightly more time depending on the size of your structure but it pays off. The keyword in the route section is opt=(qst2,calcfc). Indeed, I hardly encourage you to use the calcfc keyword in any optimization you want to run.

METHOD 2. If method 1 does not work, my next advice is to use the opt=ts keyword. For this method, the coordinates in your input file are those for the TS structure. Here is an example of the input file.

link 0
--blank line--
#p b3lyp/6-31G(d,p) opt=ts geom=connectivity freq=noraman
--blank line--
Charge Multiplicity
Coordinates of TS
--blank line--

The question that arises here is how should I get the coordinates for my TS? Well, honestly this is not a trivial task, here is where you use all the chemistry you know. For example, you can start with the coordinates of your reagents and manually get them closer. If you are forming a bond whose length is to be 1.5Å, then I suggest you to have that length in 1.6Å in your TS. Sometimes this becomes trial and error but the most accurate your TS structure is, based on your chemical knowledge, the easiest to find your TS will be. As another example, if you want to find a TS for a [1,5]-sigmatropic reaction a good TS structure will be putting the hydrogen atom that migrates in the middle point through the way. I have to insist, this method hardly depends on your imagination to elucidate a TS and on your chemistry background.

Most of the time when you use the opt=ts keyword the calculations crashes because of an error in the number of eigenvalues, you can avoid it adding noeigen to the route section; here is an example of the input file, I encourage you to use this method.

link 0
--blank line--
#p b3lyp/6-31G(d,p) opt=(ts,noeigen,calcfc) geom=connectivity freq=noraman
--blank line--
Charge Multiplicity
Coordinates of TS
--blank line--

If you have problems in the optimization steps I suggest you to ask the algorithm to calculate the force constants in every step of the optimization opt=(ts,noeigen,calcall) this is quite a harsh method because will consume long computing time but works well for small molecules and for complicated TSs to find.

Another ‘tricky’ way to get your coordinates for your TS is to run the qst2 calculation, then if it fails, take the second- or the third-step coordinates and used them as a ‘pre-optimized’ set of coordinates for this method.

By the way, here is another useful trick. If you are evaluating a group of TSs, let’s say, if you are varying a functional group among the group, focus on finding the TS for the simplest case, then use this optimized TS as a template where you add the moieties and use this this method. This works pretty well.

For this post we’ll leave it up to here and post the rest of Guillermo’s tricks and advice on finding TS structures next week when we’ll also discuss the use of IRC calculations and some considerations on energy corrections when plotting the full energy profile. In the mean time please take the time to rate, like and share this and other posts.

Thanks for reading!

 

Advertisements

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 May 26, 2016, in Computational Chemistry, Gaussian, GaussView, Models, Research, Theoretical Chemistry, Tricks, White papers and tagged , , , , , , , , , , , , , , . Bookmark the permalink. 16 Comments.

  1. Nazanin Taimoory

    Dear Dr. Barroso,

    I really appreciate for your nice post. I am wondering if its possible to send me an example of input file for reagent and product for doing opt=qst2?

    Thanks in advance,

    Best Regards,

    Nazanin Taimoory

    ________________________________

    • Well it depends on your study. What you need is to generate a structure of a reagent and thus get the corresponding coordinates for each atom, save them and then modify that molecule to make it look like the corresponding products (thereby keeping the numbering scheme which is very important), you get those coordinates too and then paste them in the corresponding section according to what Guillermo just wrote. It’s that simple!

      Have a nice day!

  2. Ashwin Gopal

    I am beginner in Gaussian and I want to find the transition state of trans to gauche conformational change of 1,3-but-2-ene . I already optimised trans and gauche using 6-31G(d,p) basis set .But everytime I run QST2 opt-freq , I get an error “Input conversion error in IntKMC”, before doing any calculations . I have already checked the labels.Can u help with this error.

    • If you already checked and are positive that the atoms labels remain consistent between both fragments, which is unlikely since this is what the error means, then try using another methodology to optimize the transition state such as opt=TS, I’m pretty sure we have already discussed it on the blog. I’m also pretty sure that your problem is related to a label mixing issue.

      I hope this helps

  3. Estimado Dr. Barroso,

    Muchas gracias tanto a usted como a Guillermo por tomarse el tiempo y dedicarse a hacer estos posts. Les envío un cordial saludo desde la Universidad Central de Venezuela. Now, in order to respect the readers I’ll present my doubt in English (:

    I’m painfuly new in this Computational Chemistry world. My thesis tutor is also quite unfamiliar with the area so we’re a little bit at lost sometimes, and we both are looking forward learning about Gaussian and Computational Chemistry so this investigation line may continue growing in our university. As a part of my graduation work I’m trying to determinate the Transition State of two Coumarin Tautomers using the QST2 method.

    The input is this:

    %NProc=4
    %Mem=350MW
    #P BHandHLYP/6-311++G(d,p)
    # opt=(qst2,calcfc) geom=connectivity freq=noraman
    # Gfinput IOP(6/7=3) Pop=full
    # Units(Ang,Deg)

    A

    0 1
    C 0 -3.007532 -0.849804 0.520509
    C 0 -3.504555 0.346457 -0.199849
    C 0 -2.037109 -1.471191 -0.183080
    C 0 -0.803541 -0.728091 0.095559
    C 0 -1.093450 0.617277 0.046711
    C 0 -2.435305 1.164188 -0.264558
    O 0 -0.148076 1.586468 0.076705
    C 0 0.575187 -1.179581 0.017585
    C 0 1.579649 -0.097940 0.458080
    C 0 1.189535 1.318680 0.070258
    C 0 2.946379 -0.400318 0.071410
    N 0 4.016359 -0.650308 -0.221226
    O 0 0.927296 -2.278342 -0.280477
    O 0 1.945210 2.183972 -0.170075
    H 0 -2.733718 -0.601512 1.538952
    H 0 -4.490834 0.603289 -0.545560
    H 0 -2.193400 -1.783933 -1.207922
    H 0 -2.475339 2.189325 -0.591938
    H 0 1.527796 -0.109845 1.550065

    B

    0 1
    C 0 -3.172957 -0.995899 -0.000104
    C 0 -3.412632 0.376629 0.000322
    C 0 -1.881291 -1.463140 -0.000454
    C 0 -0.814138 -0.562775 -0.000338
    C 0 -1.071369 0.803791 0.000066
    C 0 -2.371809 1.278839 0.000499
    O 0 -0.079467 1.718398 0.000189
    C 0 0.562955 -0.969342 -0.000068
    C 0 1.547694 -0.029292 0.000066
    C 0 1.247343 1.394030 -0.000256
    C 0 2.899271 -0.454699 0.000447
    N 0 3.945023 -0.915830 0.000517
    O 0 0.791319 -2.269037 -0.000100
    O 0 2.049208 2.266453 -0.000730
    H 0 -3.996771 -1.686735 -0.000156
    H 0 -4.425045 0.741465 0.000614
    H 0 -1.672014 -2.517194 -0.000866
    H 0 -2.538021 2.340686 0.000760
    H 0 1.729817 -2.462782 0.000075

    I believe everything’s as it should be, but when I run g09 it gives me this error:

    WANTED AN INTEGER AS INPUT.
    FOUND A STRING AS INPUT.
    B

    ?
    Error termination via Lnk1e in /opt/gaussian/g09/l101.exe at Fri Aug 19 11:40:39 2016.
    Job cpu time: 0 days 0 hours 0 minutes 0.6 seconds.
    File lengths (MBytes): RWF= 5 Int= 0 D2E= 0 Chk= 1 Scr= 1

    Now, I’ve investigated and found that it means the program is not reading “B” as the title of the second structure.

    I’ve tried changing the title, not putting the title, adding the basis set and all the other information before the second matrix, but nothing has worked and I’m quite close to tears right now.

    I’d be forever grateful if you could point me in the right direction.

    Muchísimas gracias de antemano,

    Mafer.

    • Hold the tears! This is very simple. You have to bear in mind that Gaussian is very strict about the syntax use of blank lines in order to know where something ends and something begins. Leave a single blank line before the title B, otherwise, how could Gaussian know that is a title and not a Boron atom with missing coordinates? That should help you get your calculation started, lets hope it also concludes successfully.

      Hasta pronto

      • María Fernanda Muñoz

        I’ve left the blank line but it still gives me the same error 😦 I’ve also changed the title so Gaussian doesn’t confuse the letter B with Boron, but it still gives me the same error 😥 Perhaps I should leave two blank lines instead of one? The computer with the program is in the lab so I must wait till monday to try again.

        Thank you so much for your fast and kind reply. If I manage to make it start on monday I’ll let you know… And also if it concludes successfully *crosses fingers*

        Feliz fin de semana!

    • Try the following:

      Delete the line with the B, do not replace it for any other title. After the end of the coordinates of A leave a blank line and then place the charge and multiplicity of B. I think the problem is that you don’t need to specify a title for each section of coordinates.

      Please let me know how it went but I’m pretty sure this will solve your long standing problem.

      Todo lo mejor!

      • María Fernanda Muñoz

        Dear Dr. Barroso:

        I’d like to thank you once again for your kind help and input. I was wondering if you could answer me another question. I’m working with tautomers of coumarins, a type of organic compound my tutor is studying. The thing is, using QST2 hasn’t been able to find the T.S. My tutor studied these kind of molecules a few years ago using Gaussian 03, with the QST2 protocol and everything went fine! But now, using G09, the QST2 hasn’t been able to locate the TS.

        I’ve been told my molecules are ‘too complicated’ to use QST2. They’re similar to this image (goo.gl/SQpfFH), but instead of a Methyl group, it has a hydroxi group. But I don’t understand why when using g03 it worked fine and now with g09 I’m having so many issues.

        I’m enormously thankful to you and your kind words and support.

        Feliz semana!

        María Fernanda.

      • Hello, María Fernanda

        I’m glad you’ve found the help you need in this little blog of mine.
        algorithms certainly change from one version to another and sometimes that leads to unexpected differences that aren’t obvious at first glance. In this case I could suggest you to try using G03 if you still have it available or just try give G09 a little help by using QST3 instead, this implies you have to make a very good guess of what the TS looks like but if you already have the TS structures from G03 in similar molecules you can use them as templates for your new ones.
        The structure in the link doesn’t look too complicated to me but maybe the proton transference is if the atoms are too far ahead. Is a bi-molecular mechanism possible?
        I hope this helps.

  4. cambiale el nombre de A a reactivos y B a productos, tal vez sea cuestion de sintaxis solamente.

  5. here i made a TS input file and calc run sucessfully see the below file
    # opt=qst2 b3lyp/6-31g geom=connectivity

    Title Card Required

    0 1
    C -3.65596400 -1.05789400 1.18374500
    C -2.66396500 -0.14799300 0.82775000
    C -1.78935900 -0.43189000 -0.22996300
    C -1.93139700 -1.63567300 -0.92710400
    C -2.93340600 -2.54433500 -0.57786400
    C -3.79525800 -2.25801300 0.47964400
    H -4.32598400 -0.83042300 2.00763900
    H -2.56058800 0.78983200 1.36514300
    H -1.25444900 -1.86307200 -1.74676700
    H -3.03551900 -3.47402600 -1.12975300
    H -4.57255700 -2.96418700 0.75619000
    C -0.70582400 0.54479300 -0.63638300
    C 0.81439600 0.08617100 1.33272400
    H 0.09951000 -0.65392500 1.70149200
    H 1.17543800 0.67383900 2.17982600
    C 1.97161000 -0.59904400 0.62891600
    C 1.89326200 -1.94724600 0.26259100
    C 3.14059900 0.11851400 0.33385200
    C 2.96039600 -2.57147600 -0.38763200
    H 0.99265000 -2.51229800 0.48990500
    C 4.20480900 -0.50327300 -0.31691800
    H 3.20848000 1.16287800 0.62493500
    C 4.11755900 -1.85005400 -0.67960500
    H 2.88651100 -3.61943600 -0.66387800
    H 5.10679300 0.06106200 -0.53580700
    H 4.94923400 -2.33375100 -1.18364500
    C -0.28820500 3.18451500 -0.55044600
    H 0.46815000 3.60162100 -1.22146600
    H -0.82656100 4.01109500 -0.08240000
    C 0.42100100 2.37019100 0.52800300
    O 0.68122346 3.11201875 1.51062522
    N 0.09103900 1.04019800 0.47949000
    S -1.44237800 2.07095800 -1.42423500
    H -0.05214200 0.06087600 -1.37173300
    H 0.79131015 1.62773983 0.07400891
    O 1.63727294 1.83559883 -0.00095237
    H 2.25040343 2.55191357 -0.18140787

    1 2 1.5 6 1.5 7 1.0
    2 3 1.5 8 1.0
    3 4 1.5 12 1.0
    4 5 1.5 9 1.0
    5 6 1.5 10 1.0
    6 11 1.0
    7
    8
    9
    10
    11
    12 32 1.0 33 1.0 34 1.0
    13 14 1.0 15 1.0 16 1.0 32 1.0
    14
    15
    16 17 1.5 18 1.5
    17 19 1.5 20 1.0
    18 21 1.5 22 1.0
    19 23 1.5 24 1.0
    20
    21 23 1.5 25 1.0
    22
    23 26 1.0
    24
    25
    26
    27 28 1.0 29 1.0 30 1.0 33 1.0
    28
    29
    30 31 2.0 36 1.0
    31
    32 35 1.0
    33
    34
    35
    36 37 1.0
    37

    Title Card Required

    0 1
    C -4.06433349 -1.19067187 1.33871244
    C -3.07233449 -0.28077087 0.98271744
    C -2.19772849 -0.56466787 -0.07499556
    C -2.33976649 -1.76845087 -0.77213656
    C -3.34177549 -2.67711287 -0.42289656
    C -4.20362749 -2.39079087 0.63461144
    H -4.73435349 -0.96320087 2.16260644
    H -2.96895749 0.65705413 1.52011044
    H -1.66281849 -1.99584987 -1.59179956
    H -3.44388849 -3.60680387 -0.97478556
    H -4.98092649 -3.09696487 0.91115744
    C -1.11419349 0.41201513 -0.48141556
    C 0.40602651 -0.04660687 1.48769144
    H -0.30885949 -0.78670287 1.85645944
    H 0.76706851 0.54106113 2.33479344
    C 1.56324051 -0.73182187 0.78388344
    C 1.48489251 -2.08002387 0.41755844
    C 2.73222951 -0.01426387 0.48881944
    C 2.55202651 -2.70425387 -0.23266456
    H 0.58428051 -2.64507587 0.64487244
    C 3.79643951 -0.63605087 -0.16195056
    H 2.80011051 1.03010013 0.77990244
    C 3.70918951 -1.98283187 -0.52463756
    H 2.47814151 -3.75221387 -0.50891056
    H 4.69842351 -0.07171587 -0.38083956
    H 4.54086451 -2.46652887 -1.02867756
    C -0.69657449 3.05173713 -0.39547856
    H 0.05978051 3.46884313 -1.06649856
    H -1.23493049 3.87831713 0.07256744
    C 0.01263151 2.23741313 0.68297044
    O 0.27285397 2.97924088 1.66559266
    N -0.31733049 0.90742013 0.63445744
    S -1.85074749 1.93818013 -1.26926756
    H -0.46051149 -0.07190187 -1.21676556
    H 1.17321068 2.02355211 -0.13871174
    O 2.04564243 1.96837670 -0.15591981
    H 2.40375472 2.83508421 -0.36128212

    1 2 1.5 6 1.5 7 1.0
    2 3 1.5 8 1.0
    3 4 1.5 12 1.0
    4 5 1.5 9 1.0
    5 6 1.5 10 1.0
    6 11 1.0
    7
    8
    9
    10
    11
    12 32 1.0 33 1.0 34 1.0
    13 14 1.0 15 1.0 16 1.0 32 1.0
    14
    15
    16 17 1.5 18 1.5
    17 19 1.5 20 1.0
    18 21 1.5 22 1.0
    19 23 1.5 24 1.0
    20
    21 23 1.5 25 1.0
    22
    23 26 1.0
    24
    25
    26
    27 28 1.0 29 1.0 30 1.0 33 1.0
    28
    29
    30 31 2.0 32 1.5
    31
    32
    33
    34
    35 36 1.0
    36 37 1.0
    37
    But in output file i did not observe the negative frequency if need may i post output file?

    • Dear Fazil,

      One source of error is that your input file does not request a frequency analysis at all! Unless you are doing a separate frequency calculation.
      Include freq=noraman into your route section. Also, you can modify the optimization with opt=(QST2,noeigentest)

      I hope this helps

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

%d bloggers like this: