Rep:BT1215 CP3MD Lab

Introduction

The ability to model reactions using computational simulations is something of great interest, particularly with regard to understanding experimental mechanisms or molecular properties^[1]. Modern computational software, such as Gaussian, is able to use a series of quantum mechanical calculations to optimise a reaction system and determine a lowest energy reaction pathway, as well as the predicted energies and structures of any reactants, transition states, and products involved. Within these calculations, there are two important theoretical considerations to take into account:

1. Potential Energy Surface (PES) - The potential energy of a system given as a function of molecular geometry. By changing the molecular geometry (e.g. bond lengths during a reaction), the change in energy and transition state can be computed.

2. Computational Method - Defines the approximations made within the quantum mechanics, balancing accuracy (fewer approximations lead to more precise results) with computational cost (fewer assumptions mean more resource-intensive calculations).

Potential Energy Surface (PES)

As mentioned before, the PES is a representation of the relationship between the potential energy of a system and its molecular geometry. Given that non-linear molecules have a geometrical freedom of 3N-6 modes (or 3N-5 for linear molecules), where N = the number of atoms in the molecule, the PES therefore has the same dependency, relying on the internal coordinates of the atoms in the system. In order to model a reaction PES, the number of degrees of freedom that are varied is often reduced in order to make the calculation computationally feasible. When plotted graphically, the PES represents a landscape of high and low energy configurations that correspond to turning or stationary points on the surface, with an example plot shown in Figure 1^[2]. Low energy turning points correspond to an energetically stable molecular geometry, which are mathematically characterised by having a first derivative equal to 0 (Eq 1), and a positive second derivative (Eq 2), where $V$ is potential energy and $x$ represents an atomic or reaction coordinate:

\frac{\partial V}{\partial x} = 0

Eq 1

\frac{\partial^{2} V}{\partial x^{^{2}}} > 0

Eq 2

The potential energy well surrounding the minimum can be modelled as a Simple Harmonic Oscillator, where the lowest energy point is described by a Taylor expansion. Assuming a harmonic motion modelled by Hooke's law, force constant $k$ can be equated to the second derivative in Eq 2. Since the vibrational wavenumber, $v$ , is directly proportional to $\sqrt{k}$ , it can therefore be seen that the molecular configurations corresponding to minimum turning points will only have positive vibrations^[4].

The transition state of a reaction corresponds to a saddle point which separates two local minima in the PES, and is essentially a maximum turning point along a particular reaction coordinate of the surface. This can be modelled as a barrier separating two wells, shown in Figure 2^[2]. An activation energy is required to overcome this transition state and thus interconvert between the two wells, often denoted reactants and products. The saddle point on the PES surface can be mathematically described as being a local minimum in all directions except for one, where it is equivalent to a maximum turning point. This leads to one reaction coordinate where the second derivative of potential energy is negative. Following the same key discussion above, one negative vibration frequency will be observed in the transition state, while all other frequencies observed will be positive.

It is important to note that, while this discussion of a singular force constant holds true for one-dimentional systems, polyatomic systems require more complicated treatment, since there are multiple force constants and vibrational motions within the system that can actually couple (i.e. one vibration could affect the likelihood of another happening). The general form of these force constants is shown below in Eq 3, where i and j represent degrees of freedom (or coordinates) up to 3N-6:

\frac{\partial^{2} V}{\partial x_{i} x_{j}} = k_{i j}

Eq 3

Software such as Gaussian simplifies these complicated raw motions by first converting the internal coordinates into mass-weighted coordinates to make the coupled force constants equally weighted, before creating a large Hessian matrix of all the coupled force constants. The software then diagonalises the Hessian matrix to give the decoupled, vibrational modes of the molecule^[5]. A simplified example of the Hessian matrix that is diagonalised is shown below in Figure 3, where $K_{i j}$ represents the force constants with respect to the mass-weighted coordinates.

Nf710 (talk) 11:35, 16 April 2018 (BST) Excellent understanding of how the force constants are derived. You have clearly done lots of extra reading.

Computational Method

In order to compute the potential energy surface and obtain molecular energies, Gaussian uses quantum mechanical calculations based on the Linear Combination of Atomic Orbitals (LCAO) method. As implied by the existence of the PES, the Born-Oppenheimer approximation is used to separate the electron and nuclear timescales, meaning the nuclei are effectively static on an electron timescale and can be manipulated as such. If this approximation was not true then it would be impossible to calculate the PES, since the electronic energies would constantly be affected by nuclear motion. The LCAO method is essentially a sum of atomic orbitals, $ϕ_{i}$ , which combine with a weighting coefficient, $c_{i}$ , to give the overall molecular orbital wavefunction, $ψ_{m}$ , shown below in Eqn 4:

ψ_{m} = \sum_{i}^{N} c_{i} ϕ_{i}

Eq 4

The number of atomic orbitals used to build the LCAO is called the Basis Set, which will be discussed further below. The total energy of the molecule is then calculated as a sum of the orbital energies, using the Hamiltonian Operator to solve the Schrödinger equation, where Eqn 5 shows the general form, and Eqn 6 shows the same equation but with the sum of atomic orbitals:

< ψ | \hat{H} | ψ > = E

Eq 5

\sum_{i}^{N} \sum_{j}^{N} < ϕ_{i} | \hat{H} | ϕ_{j} > c_{i} c_{j} = E

Eq 6

Gaussian represents these equations in matrix form and solves to give the molecular orbital energies, as well as the overall molecule energy. Both the method of calculation and the size of the basis set used to construct the molecular orbital are fundamental factors in determining the accuracy of the computed MOs and energy output, but also in the amount of computational power required. A variety of methods and basis sets can be used to perform optimisations and energy calculations, depending on the computational method chosen. The two methods chosen for all of the following calculations are PM6 and B3LYP (with a 6-31G (d) basis set).

PM6 is a Semi-Empirical method based on the Hartree-Fock method, however it utilises experimental data and/or DFT results (where experimental data is not available) to help speed up the calculation time. Naturally, these assumptions require the system to be modelled by the experimental data used, which may not be an accurate representation. As with the Hartree-Fock method, it also treats electrons as largely independent, and does not account for electron correlation. As a result it represents a rough approximation of the system, although its speed does make it useful for very large systems which would require far too much computational power with a more accurate method. ^[6]

B3LYP is a hybrid method based on both the Hartree-Fock and Density Functional Theory (DFT) techniques. The Hartree-Fock method is employed in the calculation of the exchange-correlation energy, while DFT is used elsewhere due to its improved efficiency. This is because it only depends on a 3-coordinate system describing the electron, and thus scales 3-dimensionally with the number of basis functions, versus the four-dimensional scaling of the Hartree-Fock method. The smaller scaling results in faster computation, though it is still a lot slower than Semi-Empirical methods. A 6-31G (d) Pople basis set used for the B3LYP calculations, where the numbers represent the basis functions used in the calculation. In general terms, a larger basis set corresponds to a more accurate calculation^[4].

To localise the transition state in the reaction systems below, three main methods are possible in Gaussian. The first method is adequate when there is previous knowledge surrounding the reaction mechanism, meaning it is possible to guess the transition state structure and then optimise it. Unfortunately this is very difficult to correctly guess, often missing the transition state. The second method has a higher level of accuracy, and involves optimising the reactants, before setting and freezing the distance between the atoms involved in the reaction. Optimisation to a minima is then performed, followed by a transition state optimisation. This approach provides Gaussian with more information surrounding the reaction pathway, and thus helps guide the optimisation to the correct result. This method was used for Exercise 1.

For Exercises 2, 3, and 4, a third, more complex method was chosen. This involves optimising the product, before breaking the bonds that are formed in the reaction. The distance between the atoms that form these bonds is set and frozen at a specific distance that resembles the transition state. The same calculations as in the second method are then used to identify and optimise the transition state structure. This method is the longest, but is also the most reliable, which was the main reason for its use in more complicated reaction systems.

Nf710 (talk) This section was excellent. Clearly lots of further reading and you backed up your discussions with some equations.

Exercise 1: Diels-Alder Reaction of Butadiene + Ethylene - PM6 Level

(Fv611 (talk) very good job overall, just a little hiccup on the bond length discussion.)

Reaction Overview

Butadiene and ethylene can react to form cyclohexene, shown below in Scheme 1. The reaction occurs via a type of [4+2] cycloaddition, commonly known as a Diels-Alder reaction, through a concerted syn addition^[7]. In the case of butadiene and ethylene, since butadiene (the diene) is more electron rich than ethylene (the dienophile), the reaction represents a normal electron demand Diels-Alder cycloaddition.

**Scheme 1 - Diels-Alder reaction between butadiene and ethylene**

Molecular Orbital Discussion

The frontier molecular orbital (MO) diagram for the reaction between butadiene and ethylene is shown below in Figure 1, and was generated using the calculated Gaussian energies of the transition state and reactants. The transition state and reactants were used in the MO diagram for clarity, since conformational changes and significant orbital mixing in the product make the connection more complicated. As a result, the energies computed for the overlapping MOs are significantly higher than those of the optimised product, especially since the transition state represents a strained conformation. It is also important to note that these energies can only be considered relative to each other due to the inaccuracy of the PM6 optimisation method.

**Figure 4 - MO diagram showing frontier orbital interactions in the Diels-Alder reaction between butadiene and ethylene.**

Interactive JMOLs of the two starting materials, as well as the orbitals involved in the bonding interaction can be seen below in Table 1, labelled with both their number, computed energy and occupancy. It is possible to toggle through the generated frontier MOs using the drop down box. Please note that they may not appear to work correctly until after all of the Jmols on the page have loaded.

Table 1 - Frontier orbitals in the Diels-Alder Reaction between butadiene + ethylene

Reactant Orbitals

Transition State Orbitals

Product Orbitals

Orbitals can form stabilising overlap interactions when their associated overlap integral is non-zero. The integral is a quantitative representation of the spatial overlap of the orbitals, and thus an integral equal to 0 is equivalent to no spatial overlap. In order for the overlap integral to be non-zero, the overall integral of the interacting orbital wavefunctions must be a symmetric function. As a result, only orbitals which are both symmetric (S × S = S) or both antisymmetric (AS × AS = S) result in a non-zero overlap integral that allows for a stabilising orbital overlap. Interaction between a symmetric and antisymmetric orbital would lead to an overall antisymmetric function (S × AS = AS) and an integral equal to 0 (representing no interaction). ^[4] This can also be seen in the MO diagram above, where only orbitals of the same symmetry interact in the Diels-Alder reaction.

Bond Length Discussion

Diels-Alder reactions involve the creation of two new σ-bonds between the diene and dienophile, while three π-bonds are lost and another one is created to complete the cycloaddition. The curly-arrow mechanism for the cycloaddition between butadiene and ethylene is shown below in Figure 5, where each carbon has been numbered so as to allow for the discussion of how each bond length changes throughout the course of the reaction.

The changes in bond lengths with the reaction coordinate are shown below in Graph 1, and were extracted from the PM6 optimised IRC reaction pathway. The typical values for bond lengths between two sp³ hybridised carbon atoms, two sp² carbon atoms (both single and double bonds), and sp³ and sp² carbon atoms is shown in Table 2, where all values were obtained from the CRC Handbook of Chemistry and Physics^[8]. The calculated bond lengths for the PM6 optimised reactants, products, and transition state can be seen in Table 3. While there are several slightly distinct values for the van der waals radius of a carbon atom, the overall consensus corresponds to a value of 1.70 Å^[9].

**Table 2 - Typical C-C bond lengths**
Bond Type	Bond Length (Å)
C-C (sp³-sp³)	1.530
C-C (sp²-sp²)	1.460
C=C (sp²-sp²)	1.316
C-C (sp³-sp²)	1.503

**Table 3 - Optimised bond length values (PM6)**
	Bond Length (Å)
Bond	Reactants	Transition State	Products
C¹-C²	1.327	1.382	1.541
C²-C³	3.413*	2.115	1.540
C³-C⁴	1.335	1.380	1.501
C⁴-C⁵	1.468	1.411	1.338
C⁵-C⁶	1.335	1.380	1.501
C⁶-C¹	3.414*	2.115	1.540

(Fv611 (talk) The Van der Waals radius is indeed 1.7A, but for each carbon, meaning the two atoms can be considered to be interacting as soon as their bond length becomes smaller than 2 x 1.7 i.e. 3.4. Incidentally, this is why the "infinite" distance for two carbons is set to be 3.414.)

It is important to note for the reactant values marked with an asterisk that as C²-C³ and C⁶-C¹ are bonds formed during the reaction, the bond length can essentially be assumed as infinite, since there is no interaction between the carbon atoms. The same bonds in the transition state have a length of 2.115 Å, which is significantly larger than the van der waal radius for carbon atoms, suggesting a very weak interaction that still does not have a significant bonding character. In the products these bonds represent C-C (sp³-sp³) bonds, however the optimised value of 1.540 Å shows significant deviation to the typical bond value of 1.530 Å. Comparison with experimentally determined bond lengths of cyclohexene showed reasonable agreement with literature values^[10], highlighting the limitation that treating bonds individually from a valence bond theory perspective does not take into account complex orbital interactions that can affect bond length.

C³-C⁴ and C⁵-C⁶ represent the sp²-sp² double bonds of butadiene in the reactants, while C⁴-C⁵ represents the sp²-sp² single bond connecting them. Their deviation from typical bond length values can be explained by delocalisation through the overlapping p-orbitals, which lengthens the double bonds (due to loss of electron density in the stabilised bonding orbitals) and shortens the single bond (as it gains partial double bond character). Throughout the reaction, the two double bonds lengthen to form sp³-sp² single bonds, while the single bond shortens to become a sp²-sp² double bond in cyclohexene.

Finally C¹-C² represents the ethylene sp²-sp² double bond in the reactants, which again lengthens to form a sp³-sp³ single bond in the cyclohexene product. Once again, the values obtained all agree strongly with literature values for the cyclohexene product^[10].

From the bond length data, the transition state appears to have a structure similar to the reactants. Following on from Hammond's postulate, it can be suggested that this reaction must have an early transition state, whereby there is little difference in structure and therefore energy between the reactants and transition state, but a large difference between the transition state and products. This suggestion was confirmed by the IRC calculation, which follows the lowest energy pathway to determine the 1D PES, similar to that seen in Figure 2. The resultant plot, as well as the animated .gif file are shown below in Table 4.

**Table 4 - IRC plot and .gif file for the Diels-Alder reaction between Butadiene + Ethylene**
IRC Plot	IRC .gif File

Transition State Vibration Discussion

As mentioned before, the Diels-Alder reaction between butadiene and ethylene has been experimentally and theoretically proven through various experimental and theoretical studies, however the concerted nature of the pathway can also be formed by looking at the vibration which corresponds to the transition state formation. As can be seen below, the vibration shows the key C¹-C⁶ and C²-C³ bonds forming at the same time, hence identifying a synchronous reaction mechanism.

Output Calculation Files

SM: File:BT1215 CIS BUTADIENE OPT PM6.LOG
IRC: File:BT1215 COMB REACT IRC PM6.LOG
PRETS OPT: File:BT1215 COMB REACT PM6 PRETS OPT.LOG
TS: File:BT1215 BEDIELSALDER ORB COMB REACT TS PM6.LOG
SM: File:BT1215 ETHYLENE OPT PM6.LOG
PROD: File:BT1215 PROD CYCLOHEXENE OPT PM6.LOG
SINGLE POINT ENERGY REACT: File:BT1215 SPE REACT PM6.LOG