I aim for two things with this post – on the one hand to produce a manual for other people to use, but also to verify what I obtained via searching, asking a more experienced user as well as receiving a reply from Gaussian support.
I will assume that the user is familiar with the general basic idea of the format of a Gaussian input file.
The objective is to find a transition state for an electrophilic aromatic substitution reaction in decane as a solvent.
The ideal initial method to obtain a transition state is a QST2 calculation which requires a geometry for the reactants as well as the products. If it fails, a QST3 calculation which also provides a guess for the geometry of the transition state may also be attempted or as a last resort, a transition state may be manually built and optimized.
Ideally the user should have optimized geometries for the reactants as well as the product(s) in a reaction, however it is possible to build the geometry “on the spot”.
In the example I provide the input geometries for m-toludine and acetone: In this case the oxygen atom has been shifted out of the way – however the hydrogen has not been (should be though) tweaked on the m-toluidine. If bonds are broken, it is advisable to lengthen the bonds slightly as for example singe bonds will be longer than double bonds. In the product, the initial molecules have been “rebounded” for the product. “Imaginary bonds”, bonds from the product in the reactants or bonds from the reactants in the product are shows as dashed lines).
%nprocshared=4
%rwf=m-toluidine-acetone-2.rwf
%nosave
%mem=750MB
%chk=m-toluidine-acetone-2.chk
#p opt=(calcall,qst2) b3lyp/6-31++g(d,p) scrf=(solvent=n-decane) nosymm
geom=connectivity
Title Card Required
0 1
C -0.66738150 0.29867781 -0.01472477
C 0.72777850 0.29867781 -0.01472477
C 1.42531650 1.50642881 -0.01472477
C 0.72766250 2.71493781 -0.01592377
C -0.66716250 2.71485981 -0.01640277
C -1.36476350 1.50665381 -0.01540677
H -1.10194426 -0.62306930 0.00047259
H 1.27728650 -0.65383519 -0.01340977
H 1.27786250 3.66708081 -0.01598277
H -2.46436750 1.50683681 -0.01558677
N -1.40248531 3.98773009 -0.01767660
H -1.16147100 4.51273462 0.79858460
H -1.16060228 4.51160224 -0.83440824
C 2.96531624 1.50654085 -0.01383691
H 3.32245021 2.00928088 -0.88825546
H 3.32144259 2.01265639 0.85904400
H 3.32205530 0.49776307 -0.01168239
C -2.94426822 -1.38737908 -0.21119110
C -2.43095250 -2.83931123 -0.21119110
H -2.78762331 -3.34370870 0.66246165
H -2.78762737 -3.34371014 -1.08484136
H -1.36095250 -2.83932442 -0.21119359
C -4.48426822 -1.38736010 -0.21119110
H -4.84094078 -1.89192796 0.66236254
H -4.84092213 -0.37854994 -0.21099580
H -4.84094128 -1.89158925 -1.08494009
O -2.63797974 -0.83944966 0.94169704
1 2 1.5 6 1.5 7 1.0 18 0.5
2 3 1.5 8 1.0
3 4 1.5 14 1.0
4 5 1.5 9 1.0
5 6 1.5 11 1.0
6 10 1.0
7 27 0.5
8
9
10
11 12 1.0 13 1.0
12
13
14 15 1.0 16 1.0 17 1.0
15
16
17
18 19 1.0 23 1.0 27 2.0
19 20 1.0 21 1.0 22 1.0
20
21
22
23 24 1.0 25 1.0 26 1.0
24
25
26
27
Title Card Required
0 1
C -1.10949687 0.47711117 1.88331619
C 0.28566313 0.47711117 1.88331619
C 0.98320113 1.68486217 1.88331619
C 0.28554713 2.89337117 1.88211719
C -1.10927787 2.89329317 1.88163819
C -1.80687887 1.68508717 1.88263419
H -1.73208670 -0.29950945 3.46519163
H 0.83517113 -0.47540183 1.88463119
H 0.83574713 3.84551417 1.88205819
H -2.90648287 1.68527017 1.88245419
N -1.84460068 4.16616345 1.88036435
H -1.60358637 4.69116797 2.69662556
H -1.60271765 4.69003560 1.06363271
C 2.52320087 1.68497420 1.88420404
H 2.88033484 2.18771423 1.00978550
H 2.87932722 2.19108974 2.75708496
H 2.87993993 0.67619643 1.88635857
C -2.60557153 -0.76810572 1.76628847
C -2.09225582 -2.22003788 1.76628847
H -2.44892663 -2.72443535 2.63994122
H -2.44893069 -2.72443679 0.89263821
H -1.02225582 -2.22005106 1.76628598
C -4.14557153 -0.76808675 1.76628847
H -4.50224410 -1.27265460 2.63984211
H -4.50222545 0.24072341 1.76648378
H -4.50224460 -1.27231589 0.89253948
O -2.55774834 -0.63242367 2.73721431
1 2 1.5 6 1.5 7 0.5 18 1.0
2 3 1.5 8 1.0
3 4 1.5 14 1.0
4 5 1.5 9 1.0
5 6 1.5 11 1.0
6 10 1.0
7 27 1.0
8
9
10
11 12 1.0 13 1.0
12
13
14 15 1.0 16 1.0 17 1.0
15
16
17
18 19 1.0 23 1.0 27 1.0
19 20 1.0 21 1.0 22 1.0
20
21
22
23 24 1.0 25 1.0 26 1.0
24
25
26
27
The calculation can then be started. Once it has successfully finished, it is advisable to verify that the molecule exhibits one imaginary frequency only and that the motion of the atoms or molecules under that frequency agree with the expected reaction. If this is the case, one should follow up with an IRC calculation to follow the reaction path to verify one has indeed found the desired transition state. (The geometry below is the transition state from the calculation above, set up for an IRC calculation.)
%nprocshared=4
%rwf=m-toluidine-acetone-2.rwf
%nosave
%mem=750MB
%chk=m-toluidine-acetone-2.chk
#p irc=(maxpoints=200,recalc=5,calcfc) b3lyp/6-31++g(d,p)
scrf=(solvent=n-decane) nosymm geom=connectivity
Title Card Required
0 1
C -1.20983500 0.26428000 0.39208600
C 0.14314000 0.04829600 0.07611800
C 0.97119600 1.12420400 -0.26335000
C 0.42483200 2.41568200 -0.28932800
C -0.92232100 2.65929200 0.02889700
C -1.73598700 1.56297000 0.36139500
H -1.80370000 -0.35728900 1.39005500
H 0.56171000 -0.95292700 0.09719200
H 1.06141800 3.25561300 -0.56193000
H -2.77794400 1.73200600 0.62392300
N -1.45199400 3.94909800 -0.04588200
H -2.29620600 4.10853800 0.48772900
H -0.78926200 4.70638200 0.05450600
C 2.43332600 0.91594500 -0.58932400
H 2.68752500 1.33957100 -1.56739900
H 3.07670000 1.40280600 0.15294900
H 2.68804800 -0.14719300 -0.60548500
C -2.49137900 -1.22149000 0.43327000
C -1.66737300 -2.42548800 0.01309500
H -2.34686200 -3.28416400 0.08629100
H -1.29025700 -2.36043500 -1.01081300
H -0.84776200 -2.60357900 0.71204300
C -3.57286200 -0.83131000 -0.55877300
H -4.30815500 -1.64705100 -0.52976700
H -4.08297000 0.08110500 -0.24873300
H -3.20072400 -0.72720100 -1.58114900
O -2.83178700 -1.17441600 1.73871300
1 2 1.5 6 1.5
2 3 1.5 8 1.0
3 4 1.5 14 1.0
4 5 1.5 9 1.0
5 6 1.5 11 1.0
6 10 1.0
7
8
9
10
11 12 1.0 13 1.0
12
13
14 15 1.0 16 1.0 17 1.0
15
16
17
18 19 1.0 23 1.0 27 1.0
19 20 1.0 21 1.0 22 1.0
20
21
22
23 24 1.0 25 1.0 26 1.0
24
25
26
27
So:
- I am aware that the initial reactants guess is not great in this calculation as it dates back to me starting to work with Gaussian.
- Is there anything I am missing that could have a significant impact on the obtained transition state (and the resulting energy barrier for the reaction)?
Answer
Mapping a Reaction Pathway: A Simple Guide
There are various ways to go about determining a reaction pathway using electronic structure theory. This post is to serve as a brief guide, highlighting a few useful techniques that may make this process more efficient. I will restrict this guide to a simple reaction mentioned below.
Preliminaries
Suppose you have a reaction of interest such as A + B -> C. In terms of relative energies we could define the electronic energies of an isolated A and an isolated B (i.e. (A)+(B)) as the separated reactants (SR) energy. Furthermore, we could define the electronic energy of C as the product (prod) energy. So now we have determined two electronic energies that can be set relative to each other in order to determine the deltaE of reaction which I will call dE which is simply (E(prod)-E(SR)). However, when it comes to reaction pathways, we are also interested in the activation energy (Ea) of the reaction. This means, for our basic reaction example, that there is one other electronic energy which we must consider, which leads us to the transition state (TS). Resolving this structure will be discussed from this point on.
Technique #1 - Finding a TS Using the QST2 & QST3 Methods
The QST2 and QST3 methods allow the user to specify two minimum energy structures lying on each 'side' of a TS, and have an algorithm determine a TS which connects the two. In our example, SR and Prod are these two structures. Now, for our simple reaction example, using QST2 is unwise as our SR is simply two molecules, A and B, separated at infinite distance. This is not a meaningful geometry in terms of a QST method. Our Prod structure would work quite well, however. The QST methods work well if both of your input geometries represent coherent molecular systems. So if you were resolving a more complicated reaction pathway with one or more intermediates, you could use the QST method to determine TS structures connecting an intermediate to another intermediate or an intermediate to a product.
The difference between QST2 and QST3 is simply that the latter allows the user to input a TS guess. If you have fully characterized both the low-lying minimum energy structures which you wish to find a connecting TS structure, you may use your own intuition and guess what the TS may or should look like. This can help speed up the process of locating a TS with QST3. However, if you do not know what the TS would look like, you can revert to using QST2 and not provide a TS guess. In my experience with these methods, it is usually more efficient to provide the program with your own user-defined internal coordinates rather than allowing the program to create its own. If you are lucky enough to have this type of job complete successfully, you will want to optimize the resulting geometry and perform a frequency analysis to determine the nature of the stationary point. Since we are looking for a TS, we would have to end up with one imaginary mode. If you don't get this, you need to go back try finding another TS.
- QST2 only requires two fully characterized minimum energy structures. These must be coherent chemical systems.
- QST3 allows the user to provide a TS guess
- Consider using your own internal coordinates rather than allowing the program to create its own
- Fully characterize the resulting structure by performing an optimization and a corresponding harmonic vibrational frequency analysis.
Technique #2 - Finding a TS via Scanning
We have already stated that SR is simply (A)+(B) where A and B are separated at infinite distance. We can look at our Prod and gain some insight as to how A and B come together to form Prod. So rather than leaving our TS guess up to some program to decide, we can be smart about the way we go about finding it by using our own, well-informed intuition. Perhaps by observing SR and Prod we make a reasonable guess that the distance between some atom on molecule A and some atom on molecule B decreases as we go from SR to Prod. Looking at the fully characterized Prod, we know the final distance between these two atoms. In SR, this distance is infinite. What we do is we take our Prod and scan across this distance, partially optimizing the geometry as we go. So if the distance between the two atoms is 1 Angstrom in the Prod, we would gradually increase this distance in increments of say 0.2 Angstroms a total of, perhaps, 10 times. Our final distance would then be 3.0 Angstroms of separation. From this can, we can look at the electronic energies and, if we guessed a good coordinate to scan over, we should have a smooth energy curve. From this we can yank the geometry with the highest electronic energy from the scan and use this as our guess to the TS.
Take your guess TS from the scan and insert it into a TS optimization job. With a bit of luck, your guess will be close to a real TS. Again, perform a frequency analysis and determine the nature of the stationary point. We must end up with one imaginary mode for our structure to be a TS.
- Scanning requires a good guess as to which coordinate to scan across; you will have to specify your own internal coordinates
- The maximum energy structure from your scan is used for a TS optimization
- This manual approach maximizes user intuition
Verifying the TS - The IRC Computation
Just because we've now characterized a TS does not mean that the TS connects, in our case, SR and Prod. We must perform an intrinsic reaction coordinate (IRC) computation. Simply, we give a program our fully characterized TS and we follow the reaction path in the forward and reverse directions. Be sure to give your IRC calculation a large number of steps. IRC calculations require force constants on your TS. Since you've already fully characterized your TS via a frequency analysis, you already have the force constants. Just feed them into your IRC job, eliminating the need to recalculate them.
Once your IRC calculation has finished, you must optimize the final structures from both the forward and reverse directions. The resulting structures should have the same energy and same geometry as the structures you used to find the TS. Of course in our simple reaction case, you would only be comparing the Prod as SR is molecules A and B at infinite distance. For this, just ensure that the IRC proceeded in such a way that A and B are separating from each other across the entire IRC pathway.
- Perform an IRC on the TS you have characterized in both the forward and reverse directions
- Optimize the final structures from the IRC computation and compare to the fully characterized structures you used to initially find your TS. If they match both energetically and geometrically, then you have found a plausible TS which connects both structures.
No comments:
Post a Comment