Module 3: The Transition State

Introduction

Tutorial

In order to become familiar with the techniques which were to be used and the outputs they generated a tutorial exercise was undertaken before the main body of the computational study was done.

This tutorial involved a study of 1,5-cyclohexadiene. 1,5-cyclohexadiene can exist principly in two different forms, the anti and the gauche conformers. These two differ in the dihedral angle around the central four carbons with the anti having a bond angle of 180° with the groups anti-periplanar to one another while the gauche has a bond angle of 0° with the groups completely eclipsed. As such we would not be unreasonable to suggest that the anti form would be lower in energy

To begin with the anti form was drawn in GaussView and optimised. This was done using the Hartree Fock approximation and a 3-21G basis set. The Hartree Fock approximation models LCAO, this is the linear combination of atomic orbitals and so the molecular orbitals are built up of simple sums of the atomic orbitals. This is a widely used approximation and so is a useful place to start. The 3-21G basis set is also a widely used basis set. This is a basic example of a split valence basis set. In this the different numbers describe the number of functions which are summed to make up the approximation for a particular type of orbital, e.g. the 3 is the number of functions used to make up the approximation for the inner orbitals.

The optimisation was set up in GaussView and then run in Gaussian. This resulted in the below results:

Anti Form


Point Group = C₂
Energy = -231.69260234

In a similar fashion the same was done for the Gauche conformer. Again results are tabulated:

Gauche Form


Point Group = C₂
Energy = -231.69166702

We can compare the results of these two conformers. We initially see the two are quite similar, they both have a point group of C₂. The main, and most important difference is the energy of the two conformers. We see that the anti isomer is in fact the lower in energy this is by a value of 2.46kJ mol^-1, which is comparable to RT. This means that we will see some conversion between forms at room temperature.

What we do notice from this is that the most stable conformer is where the two C=C bonds are as far away from each other as possible. We can envisage a further more stable conformer where the two groups are even further apart. It is important that we find the lowest energy conformer as this is the ground state with respect to which we quote relevent energies such as the activation energy and the enthalpies. This is tabulated below:

Lowest Energy Gauche Form


Point Group = C₁
Energy = -231.69266120

It is possible to compare these structures generated with reference structures given^[1]. In doing so we see that the anti structure generated corresponds to "anti1" while the gauche structure corresponds to "gauche2". Finally the lowest energy structure generated is given as "gauche3". We can also compare the energies we have obtained with reference values, these are pleasing results as the two agree, indicating that the computation done has been successful. The final observation we make is that we have indeed identified the lowest energy structure.

Moving forward we can look into conformers with different point groups. As such we can model and optimise the C_i structure. In this we optimised firstly to with a 3-21G basis set and then with a more accurate 6-31G basis set, this basis set uses more functions with which to build up the orbitals being approximated and so can be seen to be more accurate. We have also changed the method used and have now moved to using DFT, density functional theory, with the B3LYP method. This is again more advanced than Hartree Fock used before. This resulted in the following output:

C_i Point Group Optimisation

Energy = -231.68029061


Point Group = C_i
Energy = -234.61112035

What we notice immediately is the much lower energy of the result at this level of theory. This is a big reduction in energy and so is a significant factor. Overall we see this reduction in energy but no real noticeable change in the geometry or structure. We see identical bond lengths, angles and other parameters. This implies that the higher level optimisation has enabled us to better calculate the energy of the structure.

A further advancement to the calculations we have performed is to take into account the frequency and vibrations of the molecule. This is done by selecting 'Frequency' as the job type in GaussView. In doing this for the C_i point group we generate the below IR Spectrum:

The crucial observation we make here is that we see no imaginary vibrations only real ones. The significance of this is that we can now infer that we have found a minimum energy structure, if we had found imaginary vibrations we would have been at a transition state structure.

Further to this we could look further into the file and do some analysis of the structure and its relative energy. In doing so we see the following lines:

We can tabulate these results and compare them to reference^[2] data. This is done below:

Comparison Of Thermochemistry Data
Quantity	Energy	Physical Meaning
Sum of Electronic and Zero-Point Energies	-234.469339	E = E_elec + ZPE
Sum of Electronic and Thermal Energies	-234.463595	E = E + E_vib + E_rot + E_trans
Sum of Electronic and Thermal Enthalpies	-234.462651	H = E + RT
Sum of Electronic and Thermal Free Energies	-234.498545	G = H - TS

We can compare these with the reference data and we find a both a good correlation and a good fit of the computational results with the literature, this leads us to have confidence in the value of the computational modelling we have done thus far.

From this we are now ready to begin to move forward and begin to model some transition state structures.

Optimising the Chair and Boat Structures

The techniques established from the tutorial section were then used to model the chair and boat transition state structures of the Cope Rearrangement. This was done by first drawing the allyl fragment CH₂CHCH₂ and optimising this with a HF approximation and a 3-21G basis set. This was the basis for the fragments which were then modelled as combining through both a chair and boat combination.

Chair Modelling

Firstly the chair transition state was guessed. This was then optimised to the transition state by selecting this from the GaussView calculation menu. The initial optimisation failed, this was due to a poor guess of the structure and so the input structure was then redrawn more carefully to look as close as possible to the predicted transition state. We also gave an additional command in the keywords box by adding Opt=NoEigen. This was to circumvent the calculation crashing if we found more than one imaginary frequency during the optimisation.

The results of all of this was the below structure with the given vibrations.

Optimisation Of Transition State


Energy = -231.61932244
Key Vibration = -817.995

The bond lengths for the bonds being made and broken can be found in looking at this optimised structure, in both cases they correspond to 2.02Å.

The chair transition structure was then optimised using a different technique which involved the freezing of the reactive bonds and atoms, known as the reaction coordinates, while minimising the rest of the molecule. This means the rest of the molecule will in theory be at its energy minima and will hopefully reduce the complexity of the calculation which needs to be performed which can save on both time and required computational power.

The minimised structure can then be optimised to a transition state structure as before however in this case we do not need to compute the force constants as we did in the original case.

Hence we perform an initial optimisation of the structures giving the below results with the reaction co-ordinate frozen.

Minimisation of Frozen Structure


Energy = -231.61518525

We can then use GaussView to generate the required Gaussian input to minimise the structure to a transition state.

Optimisation Of Transition State


Energy = -231.61932239
Key Vibration = -817.759

We can again look at the transistion state interactively, this is the imaginary vibration. We find the motion to be very similar to the original model and with identical lengths of the bonds made and broken at 2.02Å.

Both of the optimisations give key vibrations at -818cm^-1 which is in keeping with reference data for the transition which gives such a value.

Boat Optimisation

The boat transition structure was then tackled. This was done with a further method which was the QST2 study. In this the reactants and products are supplied to Gaussian and this optimises to generate the transition state between the two.

In this modelling it soon became very apparent that the numbering and geometry of the two structures needed to be accurate. In the first attempt at optimisation we were met with failure as the two structures were too far apart to generate an imaginary transition between the two.

As such we changed the bond angles in the starting material to be much closer to the product and also looked at the numbering of the two files by changing this in GaussView. When it was apparent that the two structures were numbered identically we were then able to move forward and begin the optimisation.

We proceeded with the optimisation by generating the required structures in GaussView and then selecting Job Type as Opt+Freq and choosing to optimise to a transition state. This gave us two options from which TS(QST2) was the chosen option and so we submitted the desired calculation.

When the job finished we opened the results file and found the following data:

Optimisation Of Transition State


Energy = -231.602820248
Key Vibration = -839.91

Here we see the results for the boat optimisation. The imaginary value we generate for the key vibration is more negative than the previous value though is comparable in magnitude. We also notice that this structure is noticeably higher in energy. This is not wholly unexpected due to the widespread lack of stability in a boat structure.

Following The Reaction Path

Now we had seen the two relevent transition state it made sense to follow the path of the reaction. This is something which is possible in Gaussian where the Intrinsic Reaction Coordinate (IRC) allows us to see the pathway between a transition structure and product with minimum energy. This is done by moving between the two structures along a defined pathway of a series of points.

The number of points is a very important factor, if we have too many we waste time and efficieny as well as using excessive computational power. Worse still if we use too few points we never reach a minimum energy point which is a requirement to allow us to derive anything from the study. This was first done with 50 steps, with the results given below:

Following the Reaction Path

Input Structure	Output Structure

In looking at the output file we can see this has failed to reach a minimum. We have a number of options here. Firstly we can take the final point and optimise from here. This option was explored and resulted in the below details:

Optimisation Of Transition State 2


Energy = -230.17405957

We find we do reach a minimum with this approach and so this is a valid option. A further and more accurate option would be to restart the calculation but with more points. This was also done with the number of points doubled to take into account 100 points.

Again results are tabulated:

Following the Reaction Path 2

Input Structure	Output Structure

Immediately we see a much more likely minimised structure and so we are more happy with this result. We can now take this structure and optimise this to ensure that we have indeed culminated at the lowest energy structure:

Optimising the IRC Minimum Energy Output


Energy = -230.17405959

In looking at this we see similar energies for both approaches to getting to the minimum with both having similar energies and crucially resulting in the desired structure.

Finding the Activation Energy

The final exercise was one of the most important, that of finding the actual activation energies and comparing them with references. This was key on two fronts, most basically to see whether, as we had thought we could, if we had shown the boat TS to be higher in energy. Additionally by comparison with experimental values we could see how well our modelled data matched actual experimental study. We again tabulate results to allow rapid and simple comparison to be made.

Comparison of Activation Energies
	HF/3-21G at 298.15 K	B3LYP/6-31G* at 298.15 K	B3LYP/6-31G* at 0 K	Expt. at 0 K
ΔE (Chair) (kcal mol^-1)	45.31	33.27	31.94	33.5 ± 0.5
ΔE (Boat)(kcal mol^-1)	55.67	41.58	42.01	44.7 ± 2.0

We see a very good fit between the reference data and the calculated data and so we are happy with the accuracy of the data as it has proven reproducible. This again gives support to our observations and explanations.

Additionally we now see hard evidence of the initial theory that the boat structure is higher in energy than the chair. This is due to the high energy steric interactions which occur in the boat conformer and the fact that these are reflected in the computational calculations indicates just how useful a tool computational studies can be.

Interestingly we see different values at 0K with an increase in the energy of the boat conformer at 0K while the chair energy decreases in energy. We see as such a larger energy difference between the two. This makes sense as at 0K we have no molecular motion and no averaging of any effects and so we see the two extremes of the structures and so the effects which give them different values are more prominant.

Diels-Alder Cycloaddition

The most basic Diels-Alder reaction imaginable is the [4s + 2s] cycloaddition of cis butadiene and ethylene. The basic mechanism is given below:

This study is designed to investigate a number of aspects about Diels Alder reactions. Namely we will begin to investigate the orbital interactions; these are proposed to be two pairs of HOMO and LUMO orbitals overlapping. We will then continue and model the transition state of the molecule and then MOs of the transition state. Finally we will look into the regioselectivity of the reaction and see which of the endo or exo isomers of a particular reaction are favoured.

Looking at the Orbital Interactions

It is the intention to model this computationally to study the orbital interactions in more detail. The initial step in this was to optimise both of the reactant molecules, the summary of this is tabulated below:

Optimising the Cycloaddition Reactants
cis-Butadiene

Ethylene

We could then look at the molecular orbitals of the two molecules involved to gain insight into the key interactions and overlaps which allow reaction to proceed.

cis-Butadiene Molecular Orbitals

HOMO	LUMO
Plane Reflection: Antisymmetric	Plane Reflection: Symmetric
Axis Rotation: Symmetric	Axis Rotation: Antisymmetric

Likewise we could use Gaussview to plot the Molecular orbitals of ethylene pictorially.

Ethylene Molecular Orbitals

HOMO	LUMO
Plane Reflection: Symmetric	Plane Reflection: Antisymmetric
Axis Rotation: Antisymmetric	Axis Rotation: Symmetric

It is well defined that the principal orbital interactions are those of the HOMO-LUMO pairs which overlap with each other to form the new bonds in the product.

Previous studies have shown results for the orbitals which agree with those which we observe for our computational study and so we can link the results obtained with previous experimental details

Modelling the Transition State

It was also possible to take this further and model the transition state of the diels alder reaction. This was done by first freezing the bonds formed and optimising to a minimum and then taking this forward by then optimising to a transition state a forming bond derivatives.

Optimising Diels Alder to a Minimum


Energy = -231.48586105

Optimising the Diels Alder Transition State


Energy = -231.48963220
Key Vibration = -1067.23

Invesigating the Regioselectivity

Comparing the Exo and Endo Forms
Exo

Energy = -605.71873525
Endo

Energy = -605.72132078

These results give a observation which from first principles would be unexpected. It is widely accepted that this reaction gives the exo product however when we look at the relative energies of the two forms we see the exo form is in fact almost 7kJ mol^-1 less stable than the endo form. From this we must conclude that this reaction is under kinetic control. This has direct effects on our study of transition state structures as the kinetic product is that which is formed by the lowest energy. From this we would expect that the exo transition state will be more stable than that for the endo structure.

We can look into the transition states of the two conformers to see which is in fact lower energy. This is done in a similar way to that done before by using the frozen coordinate method. As before we created the two different isomers and first froze the bonds while optimising. The checkpoint files were then opened and the Redundant Coodinate Editor was used to create coordinates at the two bonds made in the reaction to form bond derivatives.

We then ran an optimisation to the transition state with frequency analysis. The results for this process are given below: