Third Year TS and Reactivity Lab: Tutorial

Introduction

Within this lab, you will be modeling several chemical reactions and analysing them by calculating Transition State (TS) structures. Locating transition states is an important step in confirming the mechanism of a reaction - for example, is a particular reaction stepwise or concerted? Will certain reaction pathways be favoured? In the lab, you will be analysing the geometry and the MOs of the TS structure and calculating reaction barriers to gain chemical insight into reaction mechanisms.

A simple energy profile (below, left) for a reaction is plotted using the reaction coordinate. Along the reaction coordinate, the TS (2) is the point of maximum energy, linking a set of reactants (1) and products (3). The TS is, therefore, a stationary point (i.e. a gradient of zero) with a negative second derivative along the reaction coordinate.

The energy profile is a 1-dimensional simplification of the full Potential Energy Surface (PES) of the system, where the energy is a function of only one degree of freedom - the reaction coordinate. The system has been set in equilibrium in all other degrees of freedom (minimising) and the degree of freedom relating to the reaction coordinate has been varied. In reality, a chemical system undergoing a reaction will explore more than the reaction coordinate - it can be displaced away from its equilibrium position in any degree of freedom. The potential energy can be plotted against any one or a combination of these degrees of freedom or coordinates. A 2D PES is shown on the right, resulting from varying 2 coordinates. It is an example of the dimerisation of cyclopentadiene. In Coordinate 1 (left to right), the bonding atoms are brought closer together and sp3 hybridised. In Coordinate 2 (bottom to top), the cyclopentadiene molecules are rotated about the reaction centre (twisted). This PES shows some twisting is required to prevent the atoms from getting too close to the highlighted methylene.

These coordinates are in fact the first two normal modes of the system. Recalling the number of normal modes of vibration = 3N - 6, this system with 22 atoms has 60 normal modes. You can have a look at some of the other vibrations in the example above by right-clicking and selecting a model. Models 1-20 are steps of an optimisation calculation, and models 21-80 are the 60 normal modes. You should notice that the first normal mode has a 'negative' frequency. This negative is meant to depict an imaginary number, the result of calculating the frequency of a harmonic oscillator with a negative force constant (square of a negative). At every TS there must be only one imaginary frequency as all other coordinates are minimised.

Chemical intuition plays a large part in starting from a suitable guess TS structure. Below is the general procedure for setting up and submitting Gaussian TS Optimisations. The tutorial exercises will take you through approaches for generating the TS structure to use in step 1.

TS Opt Procedure 1. Generate a guess TS structure Use one of the approaches described in the tutorial exercises to generate your guess TS structure. 2. Calculation set up In the Job Type tab: Choose: ‘Opt+Freq’. Select to optimise to 'TS (Berny)'. Set 'Force Constants' to 'Calculate at First Point' - Force constants are required for TS calculations. The algorithm must have the local curvature as an input to 'know' which direction the reaction path is. In the Method tab select the method to use: For PM6: Set the method (second dropdown box) to ‘Semi-Empirical' and choose PM6 as the method. For B3LYP/6-31G(d): Set the method (second dropdown box) to DFT and choose ‘B3LYP’. For the basis set choose 6-31 and (d) for the polarisation function. In the Additional Keywords box at the bottom, type: opt=noeigen. Press 'Submit' and follow the prompts to save the input file with an appropriate name. The calcualtion should then run, these can take some time (depending on the size of the molecule and the method being used). Only one can be run at a time, wait for it to finish or if it seems to be taking a very long time (hours!) then seek help. 3. Analysing the calculation If the job succeeds, open the log file. If it doesn't, contact a demonstrator or go to the Troubleshooting page. Hopefully, the final structure should look like the TS you were expecting. To analyse the structure: Check that it is converged, (see the Analysis section of the Appendix for a reminder). Check the vibrations. Right-click on the window and select 'Results' then 'Vibrations'. If there is more than one negative frequency then you have a higher order saddle point, if there are no negative frequencies then you have a minimum energy structure. This is likely to be because the starting structure input to the calculation was not close enough to the real TS structure on the PES. Reopen your input file and try to ammend the structure to make it closer to what the TS will resemble. Also try checking the Troubleshooting page or contact one of the demonstrators. If the structure has a single imaginary frequency then it is a TS. Check that it is the correct TS by choosing it and pressing 'Start Animation'. The vibration should correspond to the reaction coordinate (i.e. it should look like it's going between the reactants and products). IRC calculations can also be run to confirm that it is the correct TS. 4. Reoptimise with a more accurate method If the results are fine and you used PM6 then you can now reoptimise with a more accurate method if necessary. Copy the structure to a new window (Shortcut: ctrl + (C, N, V) ). Changing the method, repeat step 2 and 3.

Tutorial Exercises: Generating a guess TS

Approach 1: Straight Guess

The first way to locate a TS is to just draw a guess structure and run a straight TS Optimisation on this. This usually fails and is not recommended unless you are very sure that you are close to the TS and have a very small system.

Advantages: Fastest method

Disadvantages: Very unreliable. Requires knowledge of the TS

General Procedure 1. Draw a guess TS structure using GaussView. This can be from previously optimised fragments (with copy-and-paste). In bond-forming reactions, the distances between the reacting atoms will be somewhere in between the reactant and product distances. For example, the C-C distance in the TS of a C-C forming reaction will be somewhere in between their combined Van der Waals radii and a C-C bond length. The value of 2.2 is often used for C-C bond forming reactions. 2. Set up and run the TS Opt calculation by the procedure above

Practice: Cope Rearrangement

Cope rearrangement

The classic example of a Cope rearrangement is the [3,3]-sigmatropic rearrangement of 1,5-hexadiene. Its mechanism has been debated for a long time (stepwise/dissociative/concerted) but it is now accepted that it rearranges in a concerted, pericyclic fashion via a diradical transition state. This TS bears similarity to either the cyclohexane boat or chair conformation:

Chair Transition State		Boat Transition State

Procedure First, the allyl radical fragments (CH2CHCH2) are drawn and optimised. Open a new window and in the builder window, select Element Fragment. Select Carbon Trivalent (S-A-A) for the central atom. For the neighbouring carbons choose Carbon Trivalent (S-S-D). As a reminder, clicking on the window will place the fragment where you click, replacing dangling bonds or atoms if you click on them. The fragment can be seen in the GaussView window, with the highlighted blue atom indicating which atom will be centred (this can be changed by clicking different atoms in this window). You should now have a fragment that resembles half of a TS. In addition, at the bottom of the molecular editor window it should read 8 atoms, 23 electrons, neutral, doublet. This indicates that the fragment is a radical. Now you need to optimise this fragment. Set up the job as an optimisation job (optimise to a minimum), with the PM6 method. This should complete in seconds. We will now set up the boat TS. Copy the fragment you have generated from the log file. You should see the copied fragment in the GaussView window. Open a new window and click once to place the fragment. Place the second fragment somewhere below it. Make sure your next click doesn't add a third fragment by changing the mode to Inquire in the GaussView window (question mark icon). You want to position the two fragments to resemble the boat transition state above. Use the GaussViewkeyboard commands to help you move the two individual fragments into the appropriate geometry (check the Basic GaussView Tutorial for more help.) alt + right mouse button: rotates a single fragment alt + shift + right mouse button: translates a single fragment. Move the fragments into the boat TS positions and set one fragment to be roughly 2.2 Å away from the other, such that the terminal carbons line up (the modify bond tool in Gaussian can be useful here). Try to maintain the correct symmetry. Run a TS Optimisation on the structure following the calculation procedure above. Once the calculation has run successfully, check that you have the correct TS, with one imaginary frequency at about 900i cm-1. Now, repeat for the chair conformer. To do this with the allyl fragments, rotate one of them by 180º. You should now have two transition states, but which conformers of 1,5-hexadiene do they link? Open the PM6 optimised TS structures and run IRC calculations on them. Use Always when calculating force constants and make sure you're using PM6 as the method. For each IRC calculation, take the first and last geometry and optimise them (a minimum, not a TS of course) with HF/3-21G. Now compare geometries and energies with those on this page.

Approach 3: Starting from Reactants/Products

Sometimes we will not know enough about the form of the TS to draw a guess structure straight away. In this case, we can start from either the reactants and products, changing bond lengths to resemble the TS.

Advantages: More reliable than previous methods. Doesn't require much knowledge of the TS.

Disadvantages: Requires additional steps. Difficult if minima are far away in geometry from TS. There can be problems if you start from the products but your TS resembles the minima.

General Procedure 1. At this point, you'll need to decide whether to use the reactants or the products. As a general rule, choose whichever has the fewest molecules, so for a dimerisation reaction choose the dimer. Draw and optimise this structure. Open the resulting log file and check it's the correct structure. Copy this structure into a new window. 2. Decide which lengths and/or angles will change during the reaction. Modify them using the appropriate editing tools. Often it's enough to just change the atoms directly involved in the reaction, but if this doesn't work, try altering neighbouring atoms to reflect the change in hybridisation they might undergo. In a classic Diels Alder reaction, for example, 3 double bonds become single bonds and one single bond becomes a double bond, so it might be worth changing these values to somewhere in between. The dienophile carbon atoms become sp3 hybridised, so the neighbouring atoms might need to be distorted a little to reflect this. 3. Follow approach 2: freeze the degrees of freedom you think are involved in the reaction coordinate and optimise, then run a TS optimisation on the resulting (unfrozen!) structure.

Practice: Xylylene-SO₂ Diels Alder Cycloaddition

This is an example of a hetero-Diels Alder reaction, where SO₂ is used as a dienophile. Without knowing the direction of approach of the SO₂ or whether the barrier is early or late, it can be difficult to use draw a guess TS directly. Therefore, we'll start from the product.

Procedure Select naphthalene from the Ring Fragments. Two carbons at the end will be replaced to form the SO2 part. Replace the atoms with appropriate O and S fragments, and don't forget to add the additional hydrogen atoms with Add Valence to form xylylene. Now use opt+freq to optimise to a minimum. You'll probably notice some negative frequencies due to symmetry. In the Vibrations window, use manual displacement to break symmetry and reoptimise. In future, if you expect issues with symmetry you can try to break it before the initial optimisation to save time. Delete the bonds between xylylene and SO2 and separate them by about 2.0 Å for the C-O pair and 2.4 Å for the C-S pair. Now use the Redundant Coordinate Editor to freeze the pairs and optimise to a minimum followed by a TS calculation with opt=noeigen. Analyse the results as normal. Extra Aside from endo/exo conformers, there is an entirely different way for SO2 and xylylene to react known a cheletropic or chelotropic reaction. Try using the above to find the TS for a cheletropic reaction between the sulfur of SO2 and xylylene. This time, start with an iso-indene fragment. Compare the reaction barriers and enthalpies of reaction between the two reactions.

Appendix

Analysing Calculations

Convergence

Check that the job has properly converged. For each step in optimisation, the gradient will be calculated to determine if the system is at a stationary point. To prove that a stationary point has been found, open the log file and scroll from the bottom to find a section that looks like this:

         Item               Value     Threshold  Converged?
 Maximum Force            0.000098     0.000450     YES
 RMS     Force            0.000010     0.000300     YES
 Maximum Displacement     0.001379     0.001800     YES
 RMS     Displacement     0.000379     0.001200     YES
 Predicted change in Energy=-4.942954D-08
 Optimization completed.
    -- Stationary point found.

This job has landed on a stationary point. Compare to the following, which did not converge properly:

         Item               Value     Threshold  Converged?
 Maximum Force            0.000155     0.000450     YES
 RMS     Force            0.000021     0.000300     YES
 Maximum Displacement     0.002340     0.001800     NO 
 RMS     Displacement     0.000561     0.001200     YES
 Predicted change in Energy=-2.681147D-07

Frequency Calculation

If you haven't already done so, run a frequency calculation on the job (remember to hold parameters the same - method, basis set, solvation etc).

For a reactant or product then check that all of the vibrations are positive. For a TS, check that you have one negative vibration and that it corresponds to the reaction by visualising the vibration (Results, Display Vibrations). This is a very good indication that you have the correct TS structure.

Intrinsic Reaction Coordinate (IRC) Calculation

Once you have your TS, you can run an IRC to confirm that it is the correct structure. An IRC follows the minimum energy pathway to reactants and/or products on the PES. By checking the end structures given by the IRC to the products/reactants that you expect to form then you can confirm that the TS links the correct two minima.

IRC Calculation Procedure 1. Calculation set up Using your optimised TS structure set up a Gaussian calculation. (Note: Your starting structure must have only one imaginary frequency to run an IRC) In the Job Type tab: Choose: ‘IRC’. Keep Follow IRC in 'Both directions' unless it is a symmetric TS, then it's a wasted effort to perform the calculation in both directions, so this option can be changed to 'Forward only'. Change Calculate Force constants to 'Always'. Force constants are required as you'll be starting on a stationary point and the gradient alone isn't enough. Choosing Always (instead of once) isn't necessary but helps for difficult PESs. By default, Gaussian will take 10 steps down the PES, but this often isn't enough. Check the box 'Compute more points' and choose a higher number such as 200 (note that Gaussian will stop once converged, so if you're seeking convergence, any number larger than the number of steps it takes to converge is fine). In the Method tab make sure that you use the same method at which the TS has been optimised. This is because the IRC calculation (and similarly for a frequency calculation) operates on the PES and so must be calculated at the same level of theory as the underlying PES. E.g. if you want an IRC for PM6, use the PM6 TS. It is advised that any IRC is run using PM6. This is because IRCs take a long time if you use an expensive method and basis set (e.g. B3LYP/6-31G(d)) while calculating force constants. If you only have your TS structure optimised at the B3LYP/6-31G(d) level then running a TS optimisation on this structure using PM6 and then an IRC on this may be faster. 2. Check calculation IRCs can be problematic, especially for reactions where the PES is not well defined. Some common problems with IRCs are shown below and if you are having trouble getting them to run successfully ask a demonstrator. Otherwise, if the IRC has run until a stationary point has been found in both the reverse and forward direction then you can analyse the endpoints to confirm the TS. 3. Confirming the TS Take the first and the last step of the IRC. (I.e. the endpoints of both the reverse and the forward pathways). Optimise the structures to a minimum. Once complete open the log file in GaussView - do the optimised structures look like the products/reactants you expect? Compare them to your optimised reactants and products - are they the same? If so then your TS links the correct two minima on the PES. For a unimolecular reaction then optimising the endpoint structures will work well. When there are two or more reactants, an optimisation becomes very difficult as the PES is bumpy and the gradients are small. This can yield reactants/products in the IRC that are irrelevant to the reaction - don't use these for your reaction profiles.

Reaction Barriers and Reaction Energies

Barriers and reaction energies can be calculated at room temperature and 0K using information from the log file of a frequency calculation. You will need to run a frequency calculation for each of the optimised reactants, TS and products. Using Notepad or Notepad++, open each log file and search for a section named "Thermochemistry".

For 0K, the value to use is the "Sum of electronic and zero-point Energies".

For room temperature, the value to use is the "Sum of electronic and thermal Free Energies".

If there are multiple reactants/products then make sure that you use the optimised isolated structures (the assumption that they react from infinite separation) to calculate the correct energy barrier from.

Useful Tips and information

Comparing results of jobs

When comparing energies between structures, it's important to ensure that the energies are calculated in the same way. In every case, make sure the computational method, basis set, solvation method, convergence criteria and grid are the same. This means if you are using "int=grid=ultrafine" or "opt=tight" in one calculation, it must be used in all calculations for structures that are being compared.

Bond types are ignored by Gaussian

In a quantum mechanical calculation, all Gaussian uses from the geometry is the distances between the atoms. This means specifying a bond type is meaningless and mostly aesthetic - it won't remove two hydrogens if you change a single bond to a double bond. However, they are very important in a molecular mechanical calculation! Check that the number of atoms is correct in all calculations.

Job Monitoring

For jobs that are expected to take a long time or are taking longer than expected, it's worth opening them in GaussView to see how they are doing. In GaussView, choose File then Open. Locate the log file of the job you want to monitor and make sure 'Read Intermediate Geometries is selected and open the file. Now you can periodically refresh the job with File, Refresh. Check that the job looks like it's converging to the correct geometry.