Rep:Mod3:jg1208

Module 3: Physical Computational Chemistry by James Gammerman

Introduction

The characterisation of reaction mechanisms is still one of the most challenging problems in chemistry. Structural parameters of reactants and products involved in chemical reactions can be obtained using a variety of experimental techniques, but these methods do not provide much information about connecting pathways. On the other hand, for many molecular systems, molecular orbital theoretical methods yield structural data of at least semi-quantitative accuracy, and are even predictive for some systems. One area where theoretical methods have a clear advantage over experiment is in the characterisation of transition state structures. The characterisation of reaction potential surfaces is often possible because of the capability of complex computer programs to calculate energies, energy gradients and energy second derivatives with respect to nuclear coordinates

Fig.1. 3D Potential Energy Surface indicating a transition state relative to reactants and products

The recent availability of powerful high speed computers has extended the range of applicability of these methods. Since quantum mechanics can be used via ab initio, density functional theory (DFT) and semi-empirical approaches to study any nuclear configuration of molecules, reaction potential surfaces may be explored in intricate detail using the same computational tools as applied for the equilibrium structures of reactants and products. For systems containing more than three atoms there are too many degrees of freedom to allow the full surface to be mapped out discretely. Instead, Born-Oppenheimer potential energy surfaces are characterised by their critical points, minima and saddle points(see Fig.1 ), so that these surfaces can be represented in terms of molecular structures (minima) and transition structures (saddle points) usually supplemented with a quadratic representation of the surface around these points.

The reaction path approximation has been introduced to study systems with a large number of nuclear coordinates in terms of a unique coordinate: the distance along the path connecting reactants, transition states, stable intermediates (if any) and products on the hypersurface. The reaction path can be defined as a path of steepest descent linking the critical points. Transition structures are those points on the energy surface with only one negative eigenvalue of the Hessian matrix (H), for which the elements are the second derivatives of the energy with respect to the nuclear coordinates. These transition structures correspond to situations where the energy of the system is at a minimum for all but one independent direction where it is a maximum. This direction is called the transition vector and corresponds to the negative eigenvalue of the matrix H

In the following project, the transition structures on the potential energy surfaces will be characterised for two key reactions in organic chemistry involving large systems: The Cope Rearrangement and the Diels Alder reaction.

The Cope Rearrangement

Sigmatropic rearrangements constitute an important class of concerted pericyclic reactions governed by orbital symmetry. They involve a reorganisation of electrons during which a group attached by a sigma bond migrates to the other terminus of a conjugated pi-electron system, with a simultaneous shift of the pi-electrons. Sigmatropic rearrangements are described by stating the relationship between the reacting centres in the migrating fragment and the pi system . The order [i,j] specifies the number of atoms in the migrating fragment and in the pi system respectively.

[3,3]-sigmatropic rearrangements are very important reactions, and arguably the most useful kind are those that form new carbon-carbon bonds. An example is the [3,3]-sigmatropic rearrangement of 1,5-hexadienes – this is called the Cope rearrangement, developed by Arthur Cope in 1940. A general reaction scheme is outlined in figure 2.

Fig.2. General Reaction Scheme for Cope Rearrangement

A priori there are three possible pathways for these reactions:

1. A non-concerted pathway involving sigma-bond fragmentation resulting in two allyl radicals followed by recombination (Fig. 3a)

2. A concerted pericyclic pathway where making and breaking of sigma-bonds occur synchronously (Fig. 3b)

3. A non-concerted pathway involving sigma-bond formation to give cyclohexane-1,4-diyl (Fig. 3c)

	Fig. 3a,b,c: Possible reaction pathways for Cope Rearrangement		Fig. 4: The Chair Structure of 1,5-Hexadiene		Fig: 5: The Boat Structure of 1,5-Hexadiene

The sigma-bond fragmentation mechanism is discarded because the the activation energy for the dissociation of 1,5 hexadiene (56 kcal/mol) is greater than that of the Cope rearrangement (33.5 kcal/mol). We are therefore left with two options. There is a general consensus that the 1,5-hexadiene rearrangement passes through species that have higher symmetry than reactants or products (C₂), namely the chair (C_2h) and boat (C_2v) structures, shown in Figs. 4 and 5 respectively. The goal of this part of the report is to find the lowest energy reactants, transition states and products via optimisation of their geometries and then use this knowledge to determine the true mechanism of the rearrangement.

Optimising reactants and products

1,5-hexadiene is notable for having 2 dominant conformations - antiperiplanar and gauche. Within these, geometrically distinct orientations of the double bonds relative to the central carbon atom produce further conformations which must be taken into account when performing calculations related to reaction pathways. By optimising the geometry of a 1,5-hexadiene molecule each of these conformers can be identified.

In any calculation, the method used to optimise a structure is highly significant as it affects the length and accuracy of the calculation. In this computational experiment, a quantum mechanics-based ab initio method known as Hartree-Fock (HF) was used with the 3-21G basis set. Before proceeding, the term basis set should be explained, with reference to the information provided in the module 2 report of this computational laboratory course.
A basis set is a set of mathematical functions, normally representing atomic orbitals (AOs) - initially these AOs were used in the Hartree-Fock method and were represented as Slater orbitals, which corresponded to a set of functions which decayed exponentially with distance from the nuclei. Later, Frank Boys realised that the Slater orbitals could themselves be approximated as a linear combination of Gaussian orbitals. These Gaussian basis functions are much easier to compute, because they reduce the 4-centre integrals found in Slater orbitals to finite sums of 2-centre integrals and then 1-centre integrals. These are around 5 orders of magnitude quicker to solve than Slater orbitals, so Gaussian basis functions are much cheaper to use

There are a number of different types of basis set:

Minimal Basis Sets

These contain the minimum number of basis functions (i.e. AOs) needed for each individual atom. E.g. For H: 1s; For C: 1s, 2s 2px etc. An example of this is STO-3G, which represents three Gaussian functionals (known as primitives) per Slater Type Orbital

Split Valence Basis Sets

These are are larger basis sets and can be used for more complex systems - it increases the number of basis functions per atom. Examples include 3-21G and 6-31G, which have two basis functions of different sizes for each valence orbital. e.g. For the first orbital of H: 1s and 1s'. 3-21G is the simplest basis set which gives reasonable results. In this project is was used to loosely optimise molecules to a first approximation until bigger basis sets could be used.

Polarised basis sets

Split valence basis sets allow orbitals to change size, but not shape. Polarised basis sets add in orbitals with angular momentum beyond that needed to describe the ground state which removes this limitation. For example, polarised basis sets add d functions to carbon atoms, and even add p functions to hydrogen atoms! Two examples which will appear in this project are 6-31G** and 6-311G**. 6-31G* is a significant improvement on 3-21G due to the added polarisation of atoms, and is considered the best compromise of speed and accuracy. It has become one of the most popular basis sets

Hence for this reaction, 3-21G was chosen as the basis set as it can give reasonable results quickly.

Fig.6: Anti and Gauche Conformations

Optimisation found 10 different conformers, six of which have the gauche linkage, and four of which have the anti linkage. These conformations are shown as Newman Projections in Fig.6.

Theoretical analysis of the structures finds that anti-periplanar (app) conformations of the moelcule should be favoured over gauche ones due to favourable app interactions such as σ_(C-H)/σ*_(C-H). The ten conformers given in appendix 1 were found and are shown in table 1 along with their respective energies and point groups.

Table 1 - Anti & Gauche Conformers

Structure
Energy (Hartrees)	-231.68772	-231.69167	-231.69267	-231.69153	-231.68962	-231.68916
Point group	C2	C2	C1	C2	C1	C1

Structure
Energy (Hartrees)	-231.69260	-231.69254	-231.68907	-231.69097
Point group	C2	Ci	C2h	C1

The table shows that a) the conformers found match those given in appendix 1 and b) the lowest energy conformation of the 1,5-hexadiene is gauche 3, having an energy of -231.6927. Literature explains this stability as the result of an interaction between the vinyl proton and a pi-orbital.

Recent literature has cited HF as being too inadequate an approach for proper analysis of 1,5 hexadiene because unlike DFT Hartree-Fock does not take electron correlation into account. It therefore treats 1,5-hexadiene like n-butane, for which the energy difference between the anti and gauche confomers are an order of magnitude higher than the corresponding difference in 1,5-hexadiene. A better way to determine the true relative stability would be to use a so-called 'post-Hartree-Fock' ab initio approach, such as a Coupled Cluster CCSD(T) calculation, which accounts for electron correlation by constructing multi-electron wavefunctions using an exponential cluster operator. Alternatively, the small energy differences between the anti and gauche conformers could be included in molecular mechanics parameterisations to to facilitate an accurate calculation via this method Another way to gain more accurate results was chosen: Re-optimisation was performed using a DFT approach, with a B3LYP functional and a 6-31G(d) basis set. The anti2 and gauche3 conformers had the lowest energy using HF - how would they fare when electron correlation is taken into account? Summaries of the results are given in Figs 7 and 8.

Fig 7. Summary of anti2 DFT optimisation DOI:10042/to-12688

Fig 8. Summary of gauche3 DFT optimisation DOI:10042/to-12689

DFT clearly places the anti 2 conformer at a lower energy than gauche 3, which contradicts the HF results. However, it is better to compare the geometrical parameters such as bond lengths and angles, as the fact that DFT takes electron correlation into account will have a clear effect on the total energy. These results are compared to literature in Table 2 with reference to the optimised anti2 structure in Fig. 9.

Fig.9. DFT-optimised structure of anti2 conformer

Table 2: Comparison of differently optimized structures and literature values.

Parameter	HF/3G-21	B3LYP/6-31G*	Literature
Anti 2
C3-C4 /Å	1.553	1.548	1.536
C2-C3/Å	1.509	1.504	1.507
C1=C2/Å	1.316	1.334	1.341
Average (C-H)/Å	1.075	1.094	1.108
C2-C3-C4/°	111.4	112.7	111
C1=C2-C3/°	124.8	125.3	122.5
C2=C1-H/°	121.8	121.7	120.4
C3-C2-H/°	115.5	115.7	118.4
H-C3-H/°	107.7	106.7	107.1
Energy/a.u.	-231.69254	-234.61170

The calculated bond lengths are very similar to literature, and DFT generally gives values closest to literature, which is expected as it takes electron correlation into account. It also appears to generate the most stable structures.

The energies calculated can only be confirmed with literature by including some additional terms - this requires a frequency calculation, which can also be used to characterise the critical point, that is, to confirm that a minimum has been found. This is done by inspecting all the vibrational frequencies present to see if all are real and positive.

Fig.10 contains additional thermochemical data about the energies of the anti2 confomer following vibrational analysis at the B3LYP/6-31G(d) level of theory. The bottom four lines are most significant. In order they outline:

(i) the potential energy at 0K taking zero-point vibrational energy into account: E0 = Eelec + ZPE;
(ii) the energy at standard conditions including translational, vibrational and rotational contributions to the total thermal energy: E = E0 + Evib + Erot + Etrans;
(iii) the energy containing a correction term for RT which is particularly significant when looking at dissociation reactions: H = E + RT;
(iv) the energy including entropic contributions: G=H-TS

Fig.10: Thermochemical Energies

The IR spectrum of the anti2 conformer was also computed to confirm that a minimum was indeed reached during optimisation. All of the vibrational modes in the output file are positive, which proves that a minimum was indeed reached. The IR spectrum is shown in Fig 11. It shows peaks from C-H stretches (3000-3400cm^-1), C-H bends (670-1100 cm^-1) and a C=C double bond stretch at 1734cm^-1. This particular vibration is shown in Fig 12.

Fig 11: IR Spectrum of anti2 conformer after DFT optimisation

Fig 12: Animation of double bond vibration at 1734cm-1

Optimisation of the Chair and Boat structures

Thus far all the potential reactants and products for the Cope rearrangement have been minimised to an optimal energy. However, transition state optimisations are more difficult, because the calculation needs to know where the negative direction of curvature (i.e. the reaction coordinate) is. In this section the various ways of optimising a transition structure will be explored.

Four techniques were used. The first involved predicting an approximate 'guess' transition state structure for the chair conformer, then computing the force constant matrix (i.e. the Hessian referred to in the introduction) in the first step of the optimisation, and updating these during the course of the transition state optimisation (called 'TS Berny' after Berny Schlegel, one of Gaussian's developers). However, this will not give good results if the guess structure is far from the exact structure, as the curvature of the surface may be very different at points far removed from the transition structure. Therefore, a second technique was used which involved freezing the reaction coordinate using the keyword Opt=ModRedundant and minimising the rest of the molecule. After the molecule fully relaxed, the reaction coordinate was unfrozen and the transition state optimisation was restarted.

The third technique was applied to the boat structure and utilised Gaussian's facility for automatically generating a starting structure for a transition state optimisation based upon the reactants and products that the structure connects, known as the Synchronous Transit-Guided Quasi-Newton (STQN) Method. This method uses a quadratic synchronous transit to get closer to the quadratic region around the transition state and then uses a quasi-Newton algorithm to complete the optimisation Unlike other methods, the STQN does not require a guess for the transition structure, but rather generates a guess itself which is midway between the reactants and products, in terms of redundant internal coordinates. The method was requested using the QST2 option to the Opt keyword.

Finally, following optimisation of the boat and chair transition structures, a fourth technique known as the Intrinsic Reaction Coordinate (or IRC) method was used in order to allow prediction of which conformer the reaction paths from the transition structures will lead to. This creates a series of points by taking small geometrical steps in the direction of the PES area with the highest gradient, allowing the user to follow the minimum energy path from a transition structure down to its local minimum on a PES. IRC can itself be approached in a number of ways, and this will be discussed later in the report.

Chair structure

The first technique was applied. An allyl fragment was optimised using the HF/3-21G level of theory(DOI:10042/to-12550 ), and this was then copied to produce two optimised fragments which were positioned at a separation of around 2.20Å in such a way as to resemble a predicted transition state. This shown in Fig 13.

Fig.13. Optimised allyl fragments 2.20Å apart

Then an optimisation+frequency to TS (Berny) was run to calculate the force constant at HF/3-21G including an additional keyword (Opt=NoEigen ) to prevent the calculation from crashing if more than one imaginary frequency was detected during the optimisation. If more than on imaginary frequency appears, it shows that a higher order saddle point has been found, but not a transition structure between two minima

The optimised molecule is shown in Fig.14, along with a summary and animation of the process in Figs.15 and 16. The fragment separation (between carbons 1 and 9 in the log file) falls from 2.20Å to 2.02Å, and an imaginary vibrational mode is found at -818cm^-1, which represents the simultaneous breaking and forming of a σ bond occurring in a [3,3]-sigmatropic rearrangement. The summary shows that the root-mean-square of the force constants is 0.000016, below the tolerance of 0.0003 which is the criterion by which Gaussian assesses whether an optimisation has converged and is complete.

Fig.14 Optimised chair structure

Fig.15 Summary of optimisation

Fig.16 Animation of Imaginary Frequency(click here if animation does not work

Following this, the second technique was applied by freezing the distance between the terminal carbons of the allyl fragments to 2.20Å before optimising the structure to a minimum (DOI:10042/to-12551 ). Then the 'bond' input was changed to 'derivative' and a TS Berny optimisation was performed without any force constant calculations. Similar results to the first method were found, with a single imaginary vibration appearing at -818cm^-1, confirming the identity of the turning point as a transition state. Again, the allyl fragment separation dropped from 2.20Å to 2.02Å. The results are shown in Figs 17, 18 and 19.

Fig.17 Initial chair structure

Fig.18 Optimised chair structure by redundant method

Fig.19 Summary of redundant optimisation

Fig.20 Animation of imaginary vibration at -818cm-1 , separation decreases to 2.02A. Click here if animation does not work

We have therefore seen that both techniques optimise to the same transition structure, with the same frequencies and bond lengths. We can conclude that both are adequate ways of calculating transition structures, although it should be noted that for more complicated systems the force constant calculations take longer to solve, so for this case the second technique would be preferred as it does not require analysis of the Hessian matrix.

Boat Structure

For the boat structure the optimisation is somewhat different. The third technique using the QST2 method was applied by drawing the reactant and product opposite each other (Fig.21) and allowing Gaussian to interpolate between the two structures to try to find the transition state between them. To achieve this, the atoms in the product were renumbered to match the order in the reactant. Following this an opt+freq to TS (QST2) calculation was performed at the HF/3-21G(d) level of theory (DOI:10042/to-12552 ).However, the optimisation was not successful, as shown in Fig.22. This is confirmed in Fig 23, which shows that the RMS gradient graph did not reach zero. The transition structure looks more similar to a chair than boat structure - the reason for this is that during interpolation, the program merely traanslated the top the allyl fragment and did not consider rotation around the central bond.

Fig.21 Initial reactant and product structures

Fig.22 Failure of boat structure optimisation

Fig.23 Optimisation plants showing that RMS gradient failed to reach zero

So, starting from the optimised reactant and product in order to find the true transition state was not successful - we needed to modify the QST2 method. This was done by changing the central C2-C3-C4-C5 dihedral angle to 0^o and the inner C2-C3-C4 and C3-C4-C5 angles to 100^o (Fig.24). A TS (QST2) optmisation + frequency calculation was then run in the same way as before (DOI:10042/to-12553 ) and the results are shown in Figs 24-6.

Fig.24 Modified reactant and product structures

Fig.25 Successfully optimised boat structure

Fig.26 Imaginary frequency mode of boat transition state. Click here if animation does not work

A boat transition state was found successfully! The frequency part of the calculation also found there to be a single frequency mode at 840cm^-1, confirming the presence of a transition state. It has therefore been shown that the QST2 method utilising the reactants and products to interpolate a transition structure is useful to an extent, but has its limitations, particularly if the reactants and products which are 'plugged in' at the start are not close to the transition structure. For more difficult cases, Gaussian provides a QST3 option to 'opt', which optimises a transition state structure based on the reactants, products, and a user-provided guess for the geometry of the transition structure

Intrinsic Reaction Coordinate Methods

So the reactants and products have been optimised, and we have found the boat and chair transition structures that could exist between them. But it is still practically impossibe to predict which conformer of 1,5-hexadiene will form from the transition structure. The IRC method allows us to follow the reaction mechanism by monitoring the path of minimum energy from the transition structure to the local minimum on the PES in the direction of the product. It does this by generating a set of points which correspond to tiny changes in geometry in the direction of the PES where the gradient is steepest. If there are more points, these small changes are generated more accurately. An intial calculation is made in which the force constant is calculated once at the beginning. However, this does not optimize the geometry to a minimum, giving us three options as to how to proceed:

Method 1) We could take the last point on the IRC and run a normal minimisation on this. This is quick, but if you are not close enough to a local minimum, you may end up in the wrong minimium.
Method 2) The IRC could be redone, this time using a larger number of points until it reaches a minimum. This is more reliable than option 1, but if too many points are needed then it is possible to take the wrong direction after a while and end up at the wrong structure.
Method 3) The IRC could be redone, this time computing the force constants at every step. This is the most reliable option but is also expensive and is not always feasible for large systems.

Chair Structure

Initial Calculation

For the initial calculation, the reaction coordinate was computed in only for forward direction as the reaction coordinate is symmetrical. The number of points along the IRC was set to 50 from 6, and gave the results shown in Table 3

Table 3 - IRC data performed at n=25 points

Molecule at first Point	Molecule at last point	Summary	IRC Pathway	IRC Gradient	Animation of optimisation
					thumb
[|Log file for IRC at n=25]

The RMS gradient graph shows that the geometry has not been optimised as it is still above 0.001.Furthermore, the summary shows the total energy of the transition state to be -231.6875, which is less negative (i.e. less stable) than some of the conformers found in the first section. The three methods available to us were then used.

Method 1

Method 1 was then applied, with the results shown in table 4.

Table 4 - Last point optimisation by IRC

Calculation Summary	Optimised Transition Structure	Total energy curve	RMS Gradient	Animation of Optimisation
Log file for IRC last point optimisationDOI:10042/to-12685

Table 4 shows that the RMS gradient has dropped significantly to 0.00000069 - clearly a much better optimisation has been performed. This is confirmed by the horizontal line close to zero on the graph. Symmetrisation by Gaussview revealed that the structure has the C2 symmetry group. The energy is identical to that of the gauche2 conformer found in the first section, suggesting that the transition state obtained is in fact gauche 2. The gauche2 data is displayed below table 4 with its point group, energy and relative energy for comparison.

Method 2

The second method was then applied by restarting the IRC and specifying that 100 points be taken until it reaches a minimum. The output data, in the form of the final structure and optimisation summary, animation and graphs are given in Table 5. The RMS gradient graph is below 0.00001 and shows a smooth, horizontal line, indicating that the optimisation worked.

Table 5 - Optimisation by IRC using lasr n=100

Calculation Summary	Optimised Transition Structure	Total energy curve	RMS Gradient	Animation of optimisation
Log file for IRC Method 2 DOI:10042/to-12686

Method 3

Then method 3 - the one reported to be the most reliable - was applied by re-doing the IRC sepcifying that the force constants should be computed at all 50 steps of the optimisation. The results are shown in Table.6, and show a smooth convergence of the gradient to a value of 0.000016.

Table 6 - Optimisation by IRC computing the force constants after each of 50 steps

Calculation Summary	Optimised Transition Structure	Total energy curve	RMS Gradient	Animation of Optimisation
Log file for IRC Method 3 DOI:10042/to-12687

Comparison of three methods

The laboratory script gives a reported chair energy of -231.619322 Hartrees, which is in poor agreement with the calculated values, which are all -231.69166 to five decimal places. However, as mentioned above, method 1 gives an identical value to that of the Gauche 2 conformer, suggesting that the transition structure may resemble this. According to Hammond's postulate, this means that the chair transition state occurs late in the reaction pathway.

DFT Optimisation of the Chair Structure

A DFT optimisation at the B3LYP/6-31G(d) level of theory was then performed on the HF/3-21G TS Berny optimised structure. As mentioned earlier, DFT differs from ab initio methods such as HF by including electron correlation - in layman's terms this means that it takes into account the fact that electrons in a molecular system react to each other's motion and try to avoid each other. HF calculations consider that each electron only reacts to an averaged electron density - this fails to account for the instantaneous interactions of pairs of electrons with opposite spin. But DFT can do this, because rather than calculating a wavefunction, it derives the electron distribution (or electron density function) directly. A functional is a function of another function, so this explains the name. The results of this optimisation are given in Figs. 27-9, featuring an imaginary frequency at -569cm^-1, which confirms that a TS has been found.

Fig.27 Summary of DFT Optimisation log file

Fig.28 Imaginary vibration at 569cm-1 (click here for antimation)

Fig.29 Structure of optimised chair

Inspection of the log file allows us to conclude that a successful optimisation has been carried out because
1. The RMS gradient is below its threshold value of 0.003 and
2. The RMS displacement is below its threshold value of 0.0012

Since the DFT and HF methods are so different in their approaches to calculating electron distribution, it nonsensical to compare the energies obtained by both. However, it is possible to compare the geometrical parameters for the chair transition structure obtained by both and comparing these to literature, as in Table 7

Table 7 - Geometrical Parameters for the chair transition structures optimised by HF/3-21G and DFT B3LYP/6-31G9(d) levels of theory

Parameter	HF/3-21G	B3LYP/6-31G(d)
Terminal C-C Fragmental Separation (Å)	2.02	1.98
C-C Bond Length (Å)	1.38	1.41
C-C-C-C Dihedral Angle (o)	68.5	65.2

Very precise values for bond lengths and angles have been obtained - the range for bond lengths is 0.05Å, while the bond angles are within 3.2^o of each other. It is telling that the DFT-calculated C-C bond length is lower than HF, because this is closer to the true length of 1.54Å, as might be expected given that DFT accounts for electron correlation. The difference from the real value is likely due to the fact that the bond is not yet fully formed. The lone imaginary frequency at -569cm^-1 proves that a TS has been found. The difference between this value and the -818cm^-1 vibration found by HF may be due to the fact that the DFT optimisation brings the fragments closer together - that is, closer to the product form. Bearing in mind that imaginary frequencies are no longer present in final products, it might therefore be expected that the DFT optimisation produces an imaginary frequency closer to zero (i.e. more positive) than HF found.

Boat Structure

Now to use IRC on the boat structure! An initial calculation was run as for the chair structure described above. The results are shown in Table 8.

Table 8 - IRC data performed at n=25 points

Molecule at last point	Summary	IRC Energy Curve	IRC RMS Gradient Curve
DOI:10042/to-12700 |

The key information is in the RMS graphs. The RMS gradient is 0.0009, and is not close enough to zero to warrant a true optimisation. This is shown in the RMS gradient graph which does not converge to 0.

So, method 1 will be needed, as it was for the chair structure. The results are in Table 9.

Method 1

Table 9 - Last point Optimisation by IRC

Molecule at first point	Molecule at last point	Summary	IRC Energy Curve	IRC RMS Gradient Curve
DOI:10042/to-12701 |

The optimisation has clearly converged, as the RMS gradient is 0.000007 - very close to zero.

Method 2

As for the chair structure, the IRC was restarted, specifying that 100 points be taken until it reaches a minimum. The output data are given in table 10.

Table 10 - IRC Optimisation using 100 points

Molecule at last point	Final Energy	IRC Energy Curve	IRC RMS Gradient Curve
DOI:10042/to-12702 |

The RMS gradient appears to be heading for convergence after 12 steps, but moves away from the plateau before reaching a minimum value of 0.0001 - this is not significantly less than the threshold value of 0.0003, so the IRC analysis could be further improved by the application of method 3.

Method 3

As for methods 1 and 2, method 3 was applied in the same way as for the chair. The results are displayed in Table 11

Table 11 - IRC Optimisation computing the force constants after each of 50 steps

Molecule at last point	Summary	IRC Energy Curve	IRC RMS Gradient Curve
DOI:10042/to-12702 |

Comparison of Methods

Table 12 - Comparison of IRC calculations

Calculation	Electronic Energy (Hartrees)
Initial	-231.67892442
Method 1	-231.69266121
Method 2	-231.692653021
Method 3	-231.68557510

As table 12 shows, the the four calculations gave precise results with an average value of -231.69 to two decimal places.

Method 1 gave the lowest energy, and as stated earlier it also gave an excellent RMS gradient value of 0.000007, so we can conclude that this was the most successful method of the three. The electronic energy it produced (-231.69266) is in excellent with the script value for the gauche 3 conformer which is displayed in table 13. Comparison of this with the animation of the imaginary frequency produced by method 1 in Table 13 shows a striking resemblance too, implying that the boat TS will have the same structure as gauche 3

Table 13 - Comparison of Gauche 3 conformation and IRC method 1 Optimisation

Gauche 3 conformation	Method 1 IRC animation
	Click here if animation does not work

DFT Optimisation of Boat Structure

As with the chair structure, DFT was then used to re-optimise the boat transition states which had previously been optimised using HF/3-21G. The functional used was B3LYP and the basis set was 6-31G(d) DOI:10042/to-12722 The results are shhown in figs 30-32.

Fig.30 Summary of DFT Optimisation log file

Fig.31 Imaginary vibration at 530cm-1 (click here if antimation does not work)

Fig.32 Structure of optimised chair

The calculated electronic energy value is -231.54309307 - this is in very good agreement with the script value for the boat TS of -231.54309, which shows that an optimisation using HF followed by DFT is an accurate way to find this TS. We know that a TS was found because the only negative frequency found was at -530cm^-1 - this compares to a an imaginary frequency of -840cm^-1 found by HF. The reason for this difference is like for the chair structure - DFT optimisation brings the fragments closer together, i.e. closer to the products, and thus produces an imaginary frequency closer to zero, i.e. more positive than HF.

The bond lengths and angles were compared as for the chair structure, and these are given in Table 14

Table 14 - Geometrical Parameters for the boat transition structures optimised by HF/3-21G (log file) and DFT B3LYP/6-31G(d) levels of theory (log file)

Parameter	HF/3-21G	B3LYP/6-31G(d)
Terminal C-C Fragmental Separation (Å)	2.14	2.21
C-C Bond Length (Å)	1.38	1.39
C-C-C-C Dihedral Angle (o)	64.7	64.2

Based on the arguments we gave for the vibrational analysis, we would expect DFT to give a shorter value for the interfragmental C-C bond length, but it is in fact longer by 0.07Å. This may be because the two different approaches differ enough that in this particular case, the interfragmental distances give anomalous results, but it is hard to say.

Activation Energies of Chair and Boat

In frequency analysis there is also a section called 'thermochemistry'. This tells us information about corrections to parameters such as vibrational and rotational temperatures, as welll as the zero-point energy for the system at standard conditions, which is a correction to the electronic energy of the molecule to account for the effects of molecular vibration which persist even at 0 K

The thermochemical output for both structures optimised by DFT are shown in Figs. 33 & 34.

Fig 33 Thermochemical data for chair structure optimised by DFT/B3lYP/6-31G(d)

Fig.34 Thermochemical data for boat structure optimised by DFT/B3lYP/6-31G(d)

Inspection of Figs 33 and 34 reveals that the computed values for 'sum of electronic and zero-point energies' (chair: -231.414878; boat: -234.402342) and 'sum of electronic and thermal energies' (chair: -234.408949; boat: -234.396008) all agree with the script values to 3 decimal places. The first one tell us the potential energy at 0K taking zero-point vibrational energy into account: E0 = Eelec + ZPE, while the latter one tell us the energy at standard conditions including translational, vibrational and rotational contributions to the total thermal energy: E = E₀ + E_vib + E_rot + E_trans

The activation energy of a reaction can also be found by calculating the difference between the sums of the electronic and thermal energies of the transition states and product. We calculated these earlier for the anti 2 conformation, which DFT found to be the lowest energy conformation. So all we need to do is find the difference between these values and those found in Figs 33 & 34 and then convert the units from Hartrees to kcal/mol (1 hartree = 627.509 kcal/mol ). The results are shown in table 15, along with the script values for the corresponding HF calculations, plus the experimental values.

Table 15: Summary of activation energies (in kcal/mol)

	HF/3-21G (script)	HF/3-21G (script)	B3LYP/6-31G*	B3LYP/6-31G*	Expt.
	at 0 K	at 298.15 K	at 0 K	at 298.15 K	at 0 K
ΔE (Chair)	45.70	44.69	34.28	33.21	33.5 ± 0.5
ΔE (Boat)	55.60	54.76	41.96	41.32	44.7 ± 2.0

The results from DFT optimisation are within around 3.5 kcal/mol of the experimental data for the boat transition state at 0 K, so this is reasonably accurate. The results at 25 K are even more encouraging - the activation energy for the chair transition state is within 1.3 kcal/mol of the experimental value. The Hartree-Fock optimisations are less accurate however, always 10-12 kcal/mol away from the experimental data. This might be expected as the method is less accurate than DFT, as explained in earlier sections. It is worth noting that the activation energies for the boat transition state are consistently higher than for the chair transition state, implying that more thermal energy is needed for the reaction to occur via this mechanism. Therefore, when the reaction is under thermodynamic control (i.e. high temperatures) the system has enough energy to form an equilibrium and it can go through the boat transition structure, but when under kinetic control (i.e. low temperatures) the reaction goes through the chair transition state as this is quicker to reach. As we would expect, the table shows that the activation energies are lower for both transition states at higher temperatures, because then the molecules have more kinetic energy which can overcome the barrier to reaction than at 0 K, where they only have the zero-point energy available. The freqchk utility program was then used to obtain thermochemical data at 200K. The results are shown in Figs 35 and 36.

Fig.35 Thermochemical data for chair structure optimised by DFT/B3lYP/6-31G(d) at 200 K

Fig.36 Thermochemical data for boat structure optimised by DFT/B3lYP/6-31G(d) at 200 K

The key point is that the zero-point energy of the chair transition state is just under 8 kcal/mol lower than the boat, which corroborates the earlier data. Note that at 0 K this difference was over 8 kcal/mol, i.e. the difference in activation energies is going down, so the reaction pathway is more likely to proceed via the boat transition state than before.

Diels Alder Cycloaddition

Fig 37: FMO diagrams of two HOMO-LUMO interactions

The Diels-Alder reaction between a diene and a dienophile is one of the major synthetic reactions that are used for the formation of C-C. It's usefulness lies in the mild conditions and predictable regio- and stereospecificity of the process. For most purposes the thermal reaction requires the presence of an electron-withdrawing substituent on the alkene . The reaction is exemplified by the addition of ethylene (dienophile) to cis-butadiene (diene) as shown in Fig.37. The reaction is concerted and involves the HOMO or LUMO of one fragment interacting with the LUMO or HOMO of another fragment, resulting in a bonding and antibonding orbital. The feasibility of a given reaction depends largely on the HOMO-LUMO interactions according to the following rules

• If the HOMO of one reactant can interact with the LUMO of the other reactant then the reaction is allowed
• The HOMO-LUMO can only interact when there is a significant overlap density. If the orbitals have different symmetry properties then no overlap density is possible and the reaction is forbidden.

A quantum mechanical treatment of the transition state in this reaction allows us to derive the electronic distribution, and hence a molecular orbital analysis of the fragments in order to determine which FMO's take part in the reaction. The point groups of cis-butadiene and ethylene are C_2v and D_2h respectively. The σ_v plane through the ethylene double bond will be analysed to provide information about the orbital symmetries at hand, because this plane is also found in the transition structure and product.

A qualitative prediction of which orbitals will combine can be made by inspection of the molecules involved. Ethylene is a simple example as its double bond has a pi-bonding and pi*-antibonding orbital, which directly correspond to ethylene's HOMO and LUMO. Cis-butadiene consists of two double bonds, and these combine in two different ways to produce the HOMO and LUMO. This is shown in Fig 37, in which the ethylene is replaced with a better, analagous dienophile having an electron withdrawing group: maleic anhydride. We will learn more about this particular molecule later on. If the dienophile is substituted with substituents that have pi-orbitals that can interact with the new double bond that is being formed in the product, then this interaction can stabilise the regiochemistry of the reaction. In this part of the report, the nature of the transition structure of the Diels Alder reaction, both for cases with and without substitutents, will be studied.

Optimising the reactants

Table 16 - Molecular Orbital Analysis of Cis-Butadiene and Ethylene

Name of Molecular Orbital	Molecular Orbital Diagram	Symmetry
Cisbutadiene HOMO		Anti-symmetric
Cisbutadiene LUMO		Symmetric
Ethylene HOMO		Symmetric
Ethylene LUMO		Asymmetric
DOI:10042/to-12693 ;DOI:10042/to-12694

To calculate the molecular orbitals, the optimisation method chosen for both molecules was a semi-empirical (SE) AM1 approach - this is the fourth kind of approach we have come across after molecular mechanics, ab inito methods (such as HF) and density functional theory. SE became popular in the 1960s because the quantum mechanical calculations used in ab inito methods resulted in 3 and 4-centred integrals, which were quite difficult to solve quickly at the time. SE avoids these integrals by drawing on a library of integrals that was compiled by finding the best fit of some calculated entity like geomery or energy to the experimental values. The plugging of experimental values into mathematical procedure to get the best calculated values is called paramterisation. It is the mixing of theory and experiment that makes the method 'semiemprical' The results of the AM1 calculations are given in Table 16, along with the symmetries of the molecular orbitals.

The illustrations of the orbitals can be superimposed onto a literature molecular orbital diagram to see how they matched. As shown in Fig 38, there is a clear resemblance between the two for both ethylene and butadiene.

Fig.38 MO diagram of Cis-Butadiene	Fig.Cis-Butadiene energy levels
Fig.39 MO diagram of ethylene	Fig.Ethylene energy levels

To prove that the HOMO of the butadiene overlaps with the LUMO of the the ethylene (or vice versa) for the reaction to take place, we must investigate the transition state.

Transition State Analysis

Fig 40: Modification of transition structure

The transition state of this reaction is thought to have an envelope structure formed by drawing a bicyclo compound and removing the CH₂CH₂ fragment. This modification is shown in Fig. 40. The two fragments separated by the dotted line were then optimised using the AM1 semi-empirical method followed by recombination at a separation of 2.2 Å and TS (Berny) optimisation to generate the transition state. This was then modified using the IRC method 3 and repeated using the DFT approach for comparison DOI:10042/to-12742 , DOI:10042/to-12732 , DOI:10042/to-12741

Fig 41. Animation of Am1 Optimisation

Fig. 42 Animation of DFT Optimisation

Fig 43 Animation of IRC Method 3 (DFT) Optimisation) Click here if animation does not work

Parameter	AM1 (Å)	DFT (Å)	DFT Lit. (Å)
Terminal C-C fragmental separation	2.12	2.27	2.27
C-C Bond length in Ethene fragment	1.38	1.38	1.39
C-C bond length in diene fragment	1.38	1.38	1.38
DFT Optimisation Illustrations
Optimised Structure
Lowest real frequency at 147cm-1
Imaginary negative frequency at -956cm-1

Fig.44 - Bond Lengths for Reactants, transition structure and products using DFT

The literature bond lengths are taken from the diagram in Fig.44 There is only one imaginary vibrational frequency - at -956cm^-1 - so a transition state has indeed been found.

Analysis of the geometrical data show what we were expecting: All the C-C single bonds are 1.38Å long - short of the ideal length of 1.54Å - and all the C=C double bond is longer than the ideal length of 1.35Å. These discrepancies suggest that the transition structure found has bonds which are changing in form from sigma to pi (or pi to sigma), i.e. the transition structure lies somewhere between the reactant and products.

Regarding the distance between the terminal carbons of the two fragments, this is 2.12Å - higher than the ideal C-C bond length, but within the van der Waals contact radius of 3.40Å, which implies that there is overlap of electron density between the terminal carbons. Regarding bond angles, the H-C-H angles are all in the range 114-119^o, which is between the ideal angles for sp² (120^o) and sp³(109.5^o) hybridised carbons. This provides further evidence that the bonds are changing from sigma to pi or vice versa.

The lone imaginary vibration at -956cm^-1 confirms that there is a concerted formation of sigma bonds between the butadiene and ethene molecules due to the overlap of FMOs we saw in the previous subsection. The animation in the table below Figs 41-43 show the concerted formation and breaking of C-C bonds in the reaction, which is indeed how the parent Diels-Alder reaction proceeds The lowest frequency real vibration is found at 147cm^-1 - this corresponds to a twist about the partially-formed C-C bond.

Regarding comparison of the methods, DFT gave results in near-perfect agreement with the literature. The AM1 results are also reasonable, providing results that only differ significantly from literature for the interfragmental distance.

Molecular Orbitals of the Transition State

AM1 was then used to compute the HOMO and LUMO of the transition states for this reaction, and these are shown in Figs 45 and 46.

Fig.45 - Transition Structure HOMO

Fig.46 - Transition Structure LUMO

On the molecular orbital diagram there would be bonding and antibonding orbitals which are generated from the interaction of reactant orbitals with the same symmetry. The HOMO arises because both the HOMO of cis-butadiene and the LUMO of ethene are antisymmetric through the σ_v plane, while the LUMO arises because both the HOMO of ethene and the LUMO of cis-butadiene are symmetric. These orbitals combine to form the ones shown.

Reaction of cyclohexa-1,3-diene with maleic anhydride

Fig.47 - Reaction mechanism

Another example of a Diels-Alder reaction is the cycloaddition of cyclohexa-1,3-diene with maleic anhydride (Fig. 47). Its regioselectivity is of interest and computational techniques can be used to investigate whether the exo or endo product is preferred. It is thought that the endo product is the kinetic product, because although the exo product is more stable, there is a secondary orbital interaction in the endo form which stabilises the transition state which is therefore reached first during the reaction. Can we prove this computationally? Since DFT had been shown to give the most accurate results previously, this approach was used to investigate the regioselectivity.

Optimisation & Molecular Orbital Analysis

Firstly, the reactants' geometries were optimised using DFT at the B3LYP/6-31G(d) level of theory, as this is a popular basis set which takes electron correlation into account. After this, the optimised fragments were combined to form exo and endo transition states, after which these structures were optimised again to 'TS (Berny)' in order to find the final, optimised transition state.

Fig 48 - DFT optimised cyclohexadiene

Fig.49 Summary of cyclohexadiene optimisation DOI:10042/to-12745

Fig.50: Summary of maleic anhydride optimisation DOI:10042/to-12746

Fig.51: DFT optimised maleic anhydride

The HOMOs and LUMOs of the reactants were also computed, and these are shown in Figs 52-55.

Fig.52 - HOMO of Cyclohexadiene

Fig.53 - LUMO of Cyclohexadiene

Fig.54: HOMO of Maleic Anhydride

Fig.55: LUMO of Maleic Anhydride

For the cyclohexa-1,3-diene, the HOMO is symmetric and has a node in the middle, and the antibonding (between red and green lobes) relationships are clear to see. However, there are also bonding interactions, making this orbital weakly bonding overall. The LUMO is visibly symmetric, and the nodes make the orbital antibonding overall.
For the maleic anhydride, the HOMO is symmetric and largely bonding due to the overlap between orbitals of the same sign on the five-membered ring and the oxygen atoms. The LUMO, however, is antisymmetric, with nodes on the C=C and two C=O bonds causing it to have antibonding character.

It is important to note that the HOMO-LUMO gap for maleic anhydride is significantly larger than that for the cyclohexa-1,3-diene, and its orbitals are lower down due to the electronegativity of the oxygen. This is important, because it means that the maleic anhydride has a low lying LUMO which can easily accept electron density from the HOMO of the cyclohexadiene - much easier than vice versa. This is an allowed interaction because both of these molecular orbitals are antisymmetric. Furthermore, since the two MOs are similar in energy, their interaction energy ΔE_i is small and their stabilisation energy E_stab is large, hence this is a very favourable interaction.

Endo

Now armed with the optimised reactants, the transition state was searched for using a TS (Berny) optimisation, with initial 'guess' exo and endo structures set to an interfragmental distance of 2.1Å. The results are shown in Table 17 along with the computed HOMO and LUMO.

Table 17 - Endo DFT Optimisation

Calculation Summary	Optimised Structure	Animation of imaginary frequency at -447cm-1	HOMO at -0.24230 Hartrees	LUMO at -0.06772 Hartrees
		Click here if animation does not work
DOI:10042/to-12736

This optimisation can be deemed as successful because there the RMS gradient is 0.000014, well below the threshold value of 0.00003 , and the fact that only one imaginary frequency was found at -447cm^-1 proves that a tranisition state was indeed found. The animation of this frequency is shown in table 17.

Exo

Table 18 - Exo DFT Optimisation

Calculation Summary	Optimised Structure	Animation of imaginary frequency at -449cm-1	HOMO at -0.24214 Hartrees	LUMO at -0.07838 Hartrees
		Click here if animation does not work
DOI:10042/to-12737

We can conclude that the optimisation was successful as the RMS gradient was just 0.000072, and that a transition state was found successfully as only one imaginary frequency was found in the vibration analysis at -449cm^-1, which is shown in Table 18. This shows the bond making and breaking that occurs in formation of the transition state.

Comparison of Geometrical Parameters and Energies of the Endo and Exo Products

In order to determine the regioselectivity of the reaction we must examine the differences in geometrical and energy parameters for the two products.We start by comparing the bond lengths of the exo and endo products in table 19.

Table 19 - Geometrical Parameters for the exo and endo transition structures optimised using DFT B3LYP/6-31G(d)

Parameter	Exo	Endo
C-C Fragmental Separation (Å)	2.291	2.267
C-C Bond Length (Å)	1.479	1.480
C=C Bond length (previously single bond) (Å)	1.403	1.403
C-C Bond length (previously double bond) (Å)	1.391	1.391
Through Space Distance (C=C)-(C-C) (Å)	2.961	2.872

By and large, the exo and endo forms have the same bond lengths apart from the interfragmental separation, which is 0.02Å. In order to compare the parameters we need to clarify what is taking place in this concerted reaction: Electron density moves out of pi-orbitals and into sigma-orbitals to make two new sigma-bonds which replace two pi-bonds, as shown in Fig.47. The new double bond in the cyclohexadiene fragment (which used to be a C-C) is 1.40 Å long, which is only 0.04 Å less than a typical C=C double bond, while the new single bonds between the fragments are 1.40 Å in length - a full 0.14 Å shorter than the normal C-C single bond length.

It is also worth noting that the exo form has a larger through-space distance between the (C=O)-O-(C=O) part of the anhydride and the opposing bridge of the cyclohexadiene than in the endo form. This corroborates the literature explanation for the reaction's unusual preference for the less stable exo form: The exo product experiences less steric hindrance than the endo (as it is further away). This hindrance takes the form of an interaction known as the van der Waals interaction, which is repulsive over short distances . Hence, the exo product is stablised.

But despite its apparent (relative) instability, the endo product is the favoured result of the Diels Alder reaction. This implies that it is the kinetic product, which corresponds to the product with the lowest activation energy. Let us now examine whether the calculated energies corroborate this. The activation energy of a reaction is given by the difference between the sums of the electronc and thermal energies between the transition structure and product, and the results are displayed in tables 20 and 21.

Table 20 - Comparison of energies of reactants

Energy Parameter	Cyclohexa-1,3-diene Energy (Hartrees)	Maleic Anhydride Energy (Hartrees)	Total Energy (Hartrees)
Sum of electronic & thermal energies	-233.29092	-379.22857	-612.51939

Table 21 - Comparison of endo and exo TS energies

Energy Parameter	Exo TS Energy (Hartrees)	Endo TS Energy (Hartrees)
Sum of electronic & thermal energies	-612.487659	-612.481787
Activation Energy	0.0031731	0.0027603

As the data show, 0.0004128 Hartrees (=0.26 kcal/mol) more energy is needed to reach the exo TS than endo, which confirms that endo is the kinetic product of this reaction.

Fig.56 - Secondary orbital interactions for the case of butadiene (diene HOMO) with maleic anhydride (dienophile LUMO)

Fig.57 -FMO orbital interacations

Why is the kinetic endo product preferred in Diels Alder reactions? The answer lies in the unusual phenomenon of secondary orbital interactions, illustrated in Fig.56. The dotted lines show the usual primary interaction with overlap of orbitals with matching signs, but for the endo case there is also an additional bonding interaction (bold blue lines) between the lower lobes of two p-orbitals of the diene and the upper lobes of two p-orbitals of the anhydride. This interaction is called a secondary orbital interaction, and although it does not lead to a bond, it lowers the energy of the endo transition structure relative to that for the exo structure, where it is absent A similar example of this is the reaction of cyclopentadiene and maleic anhydride, and the secondary orbital interactions for this reaction are shown with dotted lines in Fig 57.

Molecular Orbital Analysis

The TS (Berny) optimisation also allows for generation of the TS molecular orbitals and their subsequent analysis. This is shown in Fig.58 for the endo TS and Fig.59 for the exo TS.

Figure 58: Molecular orbital analysis for exo product.

HOMO orbital of Exo Transition State.

LUMO orbital of Exo Tranisiton State.

HOMO-1 orbital of Exo Tranisiton State.

The HOMO orbital is visibly derived from overlap of the antisymmetric cyclohexa-1,3-diene HOMO and maleic anhydride LUMO. It is antisymmetric. The LUMO orbital is also anti-symmetric and since it is the low-lying LUMO of the maleic anhydride which accepts the electron density from the cyclohexadiene molecule, most of the antibonding orbital character is located on the anhydride molecule. The HOMO-1 orbital, however, is symmetric. It derives from the overlap of the LUMO of the diene with the HOMO of the dienophile, and shows no secondary orbital interaction, unlike the endo TS below.

Figure 59: Molecular orbital analysis for endo product.

HOMO orbital of Endo Transition State.

LUMO orbital of Endo Tranisiton State.

HOMO-1 orbital of Endo Tranisiton State.

The HOMO is antisymmetric and is clearly derived from the interaction of the HOMO of the diene and LUMO of the dienophile. Note that this is the same as for the exo TS HOMO. The LUMO is antisymmetric, and has antibonding character due to the large number of nodes in the centre of the molecule. But the HOMO-1 is the most revealing orbital. It derives from the overlap of the maleic anhydride HOMO and cyclohexa-1,3-diene LUMO, and tellingly there is a visible chance for overlap between the pi-orbital of the C=C on the cyclohexa-1,3-diene and the oxygen atoms on the anhydride (all red lobes). This is an illustration of the secondary orbital interaction discussed earlier to explain the favourable stability of the endo TS.

Further work: Hydrogen vs Deuterium abstraction - reaction of a chlorine atom with d₁-Dichloromethane

Introduction

Fig.60 - Potential energy surface illustrating kinetic isotope effect for hydrogen and deuterium

The kinetic isotope effect is important an important phenomenon in chemistry as it allows us to determine the nature of reaction mechanisms, for example whether a reaction is an E1 or E2 elimination. They arise because a covalent bond never stops vibrating - even in the lowest vibrational state. If it did this would be a violation of Heisenberg's uncertainty principle which states that the momentum and position of a molecule cannot be known at the same time. The minimum vibrational energy a bond can have is called the zero point energy (ZPE), and this is given by the equation E₀ = 1/2hν. In order to break a covalent bond, a certain amount of energy is needed to separate from their starting position - this energy is has to raise the vibration state of the bond from the ZPE to the point at which it breaks. Since the ZPE of a C-H bond is higher than that for a C-D bond, the C-H bond has a headstart in terms of energy. The energy needed to break a C-H bond is less than that needed to break a C-D bond, so reactions in which a C-H bond is broken proceed faster than those breaking C-D bonds, provided bond breaking is occurring in the rate-determining step. In terms of elimination reactions, the above is only true for E2 reactions, not E1, so if changing C-H for C-D changes the rate of elimination then the reaction must be E2, not E1. This concept is illustrated in Fig.60 Primary isotope effects are defined as isotope effects on reactions that break or form a bond involving the atom which has undergone a change in mass.

Fig.61 Hydrogen vs Deuterium Abstraction

d1-Dichloromethane can react with chlorine atoms in either of two ways: by hydrogen abstraction to produce HCl or by deuterium abstraction to produce DCl, as shown in Fig. 61. These two reactions can have the 'same' reaction coordinate, and pass through transition states and products that are geometrically identical apart from the relative locations of the hydrogen and deuterium due to the difference in ZPEs.

In this mini-project, in accordance with the cited source, ab initio 3-21G(*) calculations were used to determine the kinetic isotope effect for hydrogen atom abstraction from d₁-dichloromethane. The kinetic isotope effect was found by calculating the transition state for hydrogen abstraction from dichloromethane A (Fig.62), and then calculating the ZPEs for the two monodeutero derivatives of this transition state, B and C (Figs 63 and 64).

Fig.62 Transition State A

Fig.63 Transition State B

Fig.64 Transition State C

Since transition structures B and C contain a heavier deuterium atom, we would expect their vibrational frequencies to be lower than A in accordance with the expression
ν = (k/m)^1/2, and therefore for their ZPEs to be lower than A in accordance with the expression E₀ = 1/2hν. Between B and C, we might expect transition state B to have the lowest ZPE, because in C the C-D bond is being broken, which means that the deuterium's contribution to the overall mass of the molecule is less. By contrast, in B the deuterium isotope is a spectator atom and the C-D bond remains intact, so we would expect this to be the heavier molecule with the lower ZPE.

Procedure

Firstly, the geometry of dichlormethane was optimised using the HF/3-21G(*) basis set. A guess transition state for the hydrogen abstraction was then drawn using the optimised structure and constraining the lengths of the reacting bonds according to the values in Fig.65. A TS (Berny) optimisation under 3-21G(*) was then performed in order to find the optimised transition state. A vibrational analysis was then used to verify that the resulting structure was a true transition state, before the same procedure was carried out for the monodeutero derivatives B and C to find their frequencies and zero-point energies.

Fig.65 Bond Length Constraints for A

Results

The structures A, B and C were optimised usign 3-21G(*) and the results are shown in figs 65-67.

Fig. 66 Structure A Optimisation Data log file

Optimised Structure of A

Summary of Optimisation

Fig. 67 Structure B Optimisation Data log file

Optimised Structure of B

Summary of Optimisation

Fig. 68 Structure C Optimisation Data log file

Optimised Structure of C

Summary of Optimisation

The optimised structures were then used to make a guess transition structure which was itself optimised at the 3-21G(*) level of theory. The optimisation data for transition state A are shown in Fig. 69.

Fig. 69 TS (Berny) Optimisation of Transition State A Data log file

Optimised Transition State of A

Summary of Optimisation

Imaginary Frequency at -2266cm-1(click here if animation does not work

The imaginary frequency at -2266cm^-1 proves that a transition state has been found. The ZPE can be found from the thermochemistry data in the log file, and is displayed in Fig. 70.

Fig.70 Zero-point energy of A

The optimisation data for transition state B is shown in Fig.71.

Fig. 71 TS (Berny) Optimisation of Transition State B Data log file

Optimised Transition State of B

Summary of Optimisation

Imaginary Frequency at -2257cm-1(click here if animation does not work

The imaginary frequency at -2257cm^-1 proves that a transition state has been found. The ZPE data is shown in Fig.72.

Fig.72 Zero-point energy of B

The optimisation data for the transition state C is shown in Fig.73.

Fig. 73 TS (Berny) Optimisation of Transition State C Data log file

Optimised Transition State of C

Summary of Optimisation

Imaginary Frequency at -1660cm-1 (click here if animation does not work)

The imaginary frequency at -1660cm^-1 proves that a transition state has been found. The ZPE can be found in the thermochemistry data in fig.74.

Fig.74 Zero-point energy of C

The ZPEs for the three transition states are displayed in table 22.

Table 22: Zero-Point Energies for Transition Structures A, B and C

Transition Structure	Description	Zero-Point Energy (kcal/mol)
A	No deuterium isotope	15.47
B	Deuterium is spectator atom; C-D intact	13.22
C	C-D bond breaks	14.53

Our hypothesis has been proven right: The order of ZPE is B < C < A. The reason for this is that in B the C-D bond remains intact, so the bond order is higher and the molecule is therefore heavier, which corresponds to a low ZPE as
ν = (k/m)^1/2 and E₀ = 1/2hν. C has a higher ZPE because the C-D bond is breaking in the transition state, so the molecule becomes lighter and the ZPE increases. And A has the highest ZPE because it does not contain any deuterium atoms at all, which are double the mass of hydrogen atoms, so this is the lightest molecule, and therefore the highest ZPE.

We can now calculate the kinetic isotope effects for the molecules. It is best to compare transition states B and C with A, as B and C have the deuterium in different positions relative to the reaction centre: In C, the C-D bond is being broken in the transition state, so we would expect C to exhibit a primary kinetic isotope effect, while in B the deuterium is geminal to the reaction centre, so we would expect B to exhibit a secondary kinetic isotope effect, i.e. a smaller effect than C.

The kinetic isotope effects were calculated using the values for ZPEs converted into eV, and plugged into the following equation : K_i/K_j = exp^{-(ΔZPE_i-ΔZPE_j)/K_BT}

Where:
K_i/K_j = kinetic isotope effect
ΔZPE_i = Change in zero point energy for the ith derivative
ΔZPE_j = Change in zero point energy for transition state A
K_BT = Boltzmann's constant in eV K^-1 = 8.6173324×10−5
And the ZPE values were converted from kcal/mol to eV (1 kcal/mol = 0.043 eV)

Plugging in the values gave the kinetic isotope effects shown in Table 23.

Table 23: Kinetic Isotope Effects for Ratios of Transition State Rate Constants

Ratio of Rate Constants for Transition States	Kinetic Isotope Effect	Type of Effect
KB/KA	0.0224	Secondary
KC/KA	0.2044	Primary

The results show that transition state B exhibits a secondary (i.e. weaker) kinetic isotope effect as the deuterium atom is geminal to the reaction centre and acts as a spectator, whereas transition state C exhibits a primary (i.e. stronger) kinetic isotope effect as the deuterium is directly involved in the C-D bond breaking. We have therefore confirmed our hypothesis.

Conclusions

Location and characterisation of transition states are an important part of chemistry, particularly in the field of chemometrics and investigation of quantitative structure-activity relationships. As transition states cannot be seen by standard spectroscopic techniques, computational methods provide a unique opportunity for chemists to analyse these transition states through optimisations, and then use this information to predict the outcome of different reactions, along with many other parameters such as bond lengths. Frequency analysis can also be performed to provide information on a molecule's vibrational frequencies, molecular orbital composition and thermochemical data. Two reactions - the Cope Rearrrangement and the Diels Alder reaction - have been successfully studied. Ab inito methods have been used to show that the Cope Rearrangment proceeds through a stabilised chair transition state, while DFT has been used to show that the Diels Alder reaction between cis-butadiene and ethyene depends on frontier molecular orbital interactions, and the regioselectivity of cyclohexa-1,3-diene reacting with maleic anhydride depends on a secondary orbital overlap favouring the endo form.
Further work has also been carried out to calculate the zero-point energies of a molecules displaying kinetic isotopic effects, and these energies and effects have been found to be in agreement with hypothesis, bringing this computational chemistry course to a successful and satisfying conclusion.