Lab 5
Chemical Reactions

5.1 Overview

The idea of the chemical reaction is central to chemistry. One way to think of a chemical reaction is as a reorganization of a set of atoms from one or more stable structures into a different set of stable structures. In the previous lab we looked at how molecular structures could be determined through calculations. We have considered molecular structures as minima on the potential energy surface (PES). A chemical reaction thus can be thought of as a path on the PES from one minimum to another. In this lab we will consider how one might go about investigating this path on the PES and in particular the highest energy point along that path, the point commonly known as the transition state (TS).

We will look at how constraining the coordinates of a molecule can be used to force its structure to be something other than a minimum of the PES. Linear Transit calculations will be examined as a way to easily map out a reaction path by performing a number of geometry optimization calculations in a single run. The transition state will be defined and a few different ways of finding the transition state will be described.

Note: Make sure that you do not delete any of your calculations from this lab and that you can identify what each calculation is about. We will make use of the results of some of these calculations in the next lab.

5.2 Geometry Constraints

Figure 5.1: Freezing coordinates

It is often useful to be able to optimize some of the geometrical parameters of a molecule while keeping others at a chosen value. This can be done with the Model:Coordinates option. Each of the geometrical parameters displayed on the right hand side of the window has a small checkbox next to it. If a box is selected then the geometrical parameter that the box corresponds to will remain frozen during any geometry optimizationi (see figure 5.1). The value of a given coordinate can be edited by typing in the relevent box. A parameter can be frozen in either Cartesian or Internal coordinates by this method but it usually makes more sense to use Internal coordinates.

Exercise 5.1: Optimize the geometry of ethene with the C-C bond frozen at 2.5 Å. What C-H bond lengths and C-C-H bond angles do you get? Compare your results with what you obtain in an optimization where no bond length is frozen.

Use a DZ basis set and LDA functional for this exercises.

5.3 Reactants, Products and Transition States

Figure 5.2: The points and paths on a PES that characterize a reaction.

In the previous lab we talked about how a stable molecular structure corresponds to a minumum on the PES. A chemical reaction is the transformation of a molecule or molecules into another different molecule or molecules. To a theoretician, the reactants and products of a reaction correspond to two different minima on a PES and the reaction is a path between the two minima (See figure 5.2).

The energies of these two minima can be used to calculate the energy of the reaction

ΔE = E_prod - E_react

(5.1)

ΔE is not the same as the reaction enthalpy ΔH or the reaction free energy ΔG but, as we shall see in the next lab, it is related to these quantities.

A path from reactants to products on the PES is a curved line through multidimensional space. There are infinitely many paths from the reactants to products. We are going to focus our attention on the paths that are minimum in energy in all directions except the one direction that leads from reactants to products (see figure 5.2). This special path we shall call a minimum energy path (MEP). The MEP is a one-dimensional curve along the PES. The progress along the MEP can be characterized by a single parameter that we shall call the reaction coordinate.

The highest energy point along a MEP, the transition state (TS) is very important. Once the energy of the transition state is known it can be used to calulate the barrier that must be overcome for the reaction to take place

ΔE^f = E_TS - E_react

(5.2)

If we know the barrier of a reaction we know something about the kinetics of that reaction. The most common theory used to calculate kinetic rate constants iis transition state theory. Transition state theory says that the rate constant of a reaction can be written as

k = C e^{-[(ΔE^f)/RT]}

(5.3)

where R is the gas constant, T is the temperature and C depends on the partition functions of reactants, products and the transition state. Thus, if C and ΔE^f could be calculated then transition state theory could be used to evaluate reaction rate constants.

Further Information: Transition state theory is also known as activated complex theory. For more information see Physical Chemistry by R.S. Berry, S.A. Rice, J. Ross, Oxford University Press, http://www.rpi.edu/dept/chem-eng/Biotech-Environ/Projects00/enzkin/transition.htm, http://www.engin.umich.edu/~CRE/03chap/html/transition/ and http://en.wikipedia.org/wiki/Transition_state.

So far, the discussion of this section considers an elementary step of a reaction and only one possible mechanism. Reactions often have more than one elementary step. This corresponds to jumps between several minima on the PES with the first and last minima being the reactants and products and the other minima being reaction intermediates. Alternative mechanisms correspond to different minimum energy paths. Each of these paths has its own transition state.

5.4 The Reaction Coordinate

The reaction coordinate is much like any other coordinate we have been discussing so far. A bond length is another one dimensional property. The reaction coordinate can be described just like a geometric coordinate. It can be given a certain value at the beginning of a reaction and a second value at the end of the reaction. The values in between correspond to the reaction taking place. The reaction coordinate thus can be used to characterize the progress of the reaction (figure 5.3).

Figure 5.3: A reaction where a bond forms between atoms A and B described by an approximate reaction coordinate. In this case that coordinate is the distance from atom A to atom B.

In general, all the atoms of a molecule will move around as a reaction takes place. Often, one particular geometrical coordinate will be the most important. An obvious example of an possible important coordinate is the distance between two atoms that form a bond in the course of a reaction. In some cases we can force the reaction of interest to occur over the course of a calculation by varying this one geometrical coordinate. If we can do this, then we can approximate the true reaction coordinate by this geometrical coordinate. It is only approximate because in general when we force a reaction to occur by varying a single geometrical parameter, the reaction path that is followed will not be the MEP and the highest energy point along this reaction path will not be the transition state. However, we hope that the path followed will be similar to the MEP and the highest energy point will be similar to the TS.

Example: In the simple acid-base reaction of H₂O and BH₃ a good choice of geometrical parameter to take as the approximate reaction coordinate is the distance between the B and O atoms.
Example: The HCN molecule can undergo an isomerization reaction where the H atom moves from the C atom to the N atom to give CNH. The approximate reaction coordinate of this reaction can be chosen to be the H-C-N angle.
Example: Ethene will bind to a metal atom through its π electrons. In this case a simple atom to atom distance is not a useful choice of approximate reaction coordinate. A better choice would be the distance between the metal atom and the midpoint between the two carbon atoms of ethene.
Exercise 5.2: Suggest an approximate reaction coordinate to describe the reaction where a molecule of carbon monoxide binds to the tungsten atom of a metal complex through the lone pair of electrons on the carbon atom of CO.
Exercise 5.3: Suggest an approximate reaction coordinate to describe the conversion of Z-1,2-dichloroethene to E-1,2-dichloroethene.

5.5 Making reactions go with Geometry Constraints

We now have a simple way to describe the progress of some reactions using an approximate reaction coordinate corresponding to a suitable geometrical coordinate. From the discussion in section 5.2 we know that we can freeze any geometrical parameter at any value we choose. By choosing reasonable values of a suitable reaction coordinate we can therefore describe a molecule at various points along a reaction path. We model a reaction taking place.

Example: If we first optimize the geometry of the H₂O-BH₃ complex and then force the O-B distance to become greater and greater while optimizing all other geometrical parameters we can describe the breaking of the O-B bond.
Exercise 5.4: We are going to consider the conversion of HCN into CNH. Optimize the geometry of HCN molecule. Freeze the H-C-N angle at 150 degrees and reoptimize the geometry. Repeat the previous step at 120 and 90 degrees. Save a picture of the geometry of the molecule at the end of each geometry optimization. Record the energy of the system at the end of each geometry optimization.

Use a DZP basis set and an LDA functional for this exercise and for all later exercises concerning HCN.

5.6 Linear Transit

The method of investigating a reaction as described in the previous section can obviously get very tedious if you are interested in a large number of different values of your chosen approximate reaction coordinate or if you don't know what values of the reaction coordinate are interesting to you and you would like to survey a range of values.

A series of calculations where one geometric variable is varied and all others are optimized for each value of the frozen variable has been automated. Such a calculation is called a Linear Transit Calculation. In a linear transit calculation starting and ending values for a given bond length, bond angle or dihedral angle are specified along with how many steps are required to go from the start to the finish. If the starting value is r_s, the finishing value is r_f and the number of steps is n then n+1 geometry optimizations are carried out with the geometrical parameter increased by (r_f-r_s)/n at each new geometry optimization. Thus in the first linear transit step the geometric parameter is r_s, in the second it is r_s+(r_f-r_s)/n, in the third it is r_s+2(r_f-r_s)/n and so on until r_f.

Example: If we run a linear transit calculation on the BH₃-H₂O complex where we specify that the parameter to be varied is the B-O bond distance and that it will start at 2.0 Å finish at 4.0 Å and have 10 steps then 11 geometry optimizations will be performed with the B-O bond distance frozen at 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8 and then 4.0 Å.

Figure 5.4: Linear transit windows.

In order to run a linear transit calculation we must first go back to the Main Options window (see figure 5.4) From the Main Options window Linear Transit should be selected from the option list next to Preset (figure 5.4).

Using the Coordinates option from the Model menu ensure that you are working with internal coordinates.

Choose the Linear Transit option from the Model menu. At the top of the Linear Transit window is a box with the number 5 in it labelled "Number of LT points." This is our number of steps n and the default value is 5.

You now need to specify which geometrical parameter is to be stepped n times and what its starting and finishing values should be. This can be a little complicated so we will describe how it is done stepwise.

First decide what geometrical parameter you want to vary.
You can only vary geometrical parameters that are defined by the internal coordinates that ADFINPUT is using. To see what internal coordinates you have, choose the Coordinates option from the Model menu and make sure that you have the internal coordinates option set.
Examine the internal coordinates and try and find the geometrical parameter that you want.
The desired geometrical parameter may not be there. There is usually many possible z-matrices that can be used to describe a given molecule. It is usually possible get ADFINPUT to give you a z-matrix that includes the parameter that you want by reordering the atoms. In general, putting the important atoms together at the top of the z-matrix in the correct order will work.
Once you have found the geometrical parameter of interest in the z-matrix, determine the atom that is associated with that parameter in the z-matrix. Select that atom in the displayed structure of the molecule.
Select the Linear Transit option from the Model menu again. Included now in the linear transit information should be the line of the z-matrix corresponding to the atom that you have selected. Below the line of the z-matrix the geometrical parameters specified by this lines of the z-matrix will be listed as Related Parameters. Next to each of these related parameters are two empty boxes. Enter in these two boxes the desired values of r_s and r_f.
The approximate linear transit coordinate must be activated by pressing the Add button next to the boxes for r_s and r_f.

Once you have followed these steps you should now be ready to run a linear transit calculation.

You can visualize the progress of a linear transit with the MOVIE program in the same way that you could follow a geometry optimization. See section 4.4.1 in the last lab.

Exercise 5.5: Convert HCN to NCH by varying the HCN angle from 180 degrees to 0 degrees in 19 steps. The symmetry of the molecule changes whenever the angle changes from 180 or 0 with the point group changing from C_infv to C_s. This confuses ADF so help it out by setting the symmetry to C_s. You do this through the menu Details: Symmetry. Once this is selected look for the drop-down menu Symbol in the righthand side window. It should be set to AUTO. Tell ADF that you want C_s symmetry by choosing C(s).

5.6.1 The Reaction Energy Profile from a Linear Transit Calculation

Figure 5.5: Energy and geometries from a linear transit calculation.

The curve described by the energy of a system as a function of reaction coordinate (the reaction energy profile) is very useful. The end points give the energy of the reactants and the products. The maximum is the energy of your approximation to the transition state. If other maxima and minima are present in the reaction energy profile then the section under consideration is not a true elementary step and intermediates and other transition states have been identified.

So far we have been looking at the geometries of the molecules as a function of the approximate reaction coordinate. The energy of the system as a function of the reaction coordinate can be found in our detailed output. To find a summary of the energies of a linear transit calculation first open the OUTPUT program using SCM:OUTPUT. If necessary, locate and open the .out file corresponding to your linear transit calculation. Choose the LT Path option from the Other Properties menu. At the top of the window the energy of the system in each linear transit step are given.

The highest energy linear transit step is the best approximation to the transition state provided by the calculation.

Exercise 5.6: Take the output of the HCN linear transit calculation that you ran in the previous section and plot the energy of the reaction as a function of the reaction coordinate. At what value of your chosen approximate reaction coordinate is the highest point on the curve? What is the barrier of the reaction? What is the overall reaction energy? Give all your energies in kJ/mol and recall that 1 a.u. of energy is equal to 2625.5 kJ/mol.

5.7 A More Refined Transition State

Figure 5.6: Read the linear transit coordinates from the detailed output and type them into ADFINPUT.

It was noted on the previous section that the highest energy linear transit step is the best approximation to the transition state provided by that calculation. How good that approximation is will depend on how good the chosen approximate reaction coordinate is and how many linear transit steps there was. The transition state can always be improved by including nore linear transit steps. With more linear transit steps the value of the reaction coordinate will change less between steps. Smaller steps means that it is possible to determine more accurately the value of the approximate reaction coordinate that has the highest energy. This approach quickly becomes very expensive in terms of computer time however because you will be running a lot of calculations.

A more efficient approach is to start from the best approximation provided by a reasonable linear transit calculation and then search for the maximum in energy by hand. Numerous ways to perform this optimization exist. The simplest method is to start from the highest energy step of the linear transit calculation and vary the reaction coordinate by small amounts until a maximum in energy is found. More sophisticated approaches make use of the gradient of the energy along the reaction coordinate. The energy gradients for all geometrical parameters can be found in the detailed output through the Iterations:Energy Gradients option of the OUTPUT program.

Figure 5.7: Find the important geometry with MOVIE and bring it into ADFINPUT with the File: Update Geometry in ADFinput menu option.

From all the geometries of the linear transit calculation we are interested in the one that is closest to the transition state. How do we choose that geometry for further study? A tedious way is to find the important geometry in the detailed output and then modify the geometrical parameters displayed by the Model: Coordinates to match what you want (figure 5.6). This method could take a very long time if your molecule is large. Alternatively, you could follow the progress of the linear transit with the MOVIE program, find the frame that corresponds to the geometry that you want and import that geometry into ADFINPUT. Finding the important frame could take a little while. A linear transit calculation is essentially a series of geometry optimizations. ADFMOVIE shows all of the geometries of each optimization and not just the final geometry for each linear transit point. You will know that you are at the frame that you want because the notice "Converged Geometry" and the correct energy will appear along the bottom of the window (figure 5.7). Once you have found the right frame the command File: Update Geometry in ADFinput will take that geometry from MOVIE into ADFinput.

Exercise 5.7: Starting from the HCN linear transit calculation, find the value of the reaction coordinate that gives the highest energy. Use whatever optimization method you desire. Your final energy should be accurate to 1 degree and 0.0005 a.u. (about 1 kJ/mol).

5.8 Transition State Search

Figure 5.8: A direct transition state search.

Rather than using an approximate reaction coordinate and constraints of geometrical parameters to find an approximation to the TS it is possible to get ADF to simply directly search for the real TS.

A direct TS search is achieved by choosing the Transition State Search option from Preset in the Main Options window (see figure 5.8). Some further options for the TS search can be modified through the Details:Task:Transition State option once a search has been chosen. Once a TS search has been run the progress of the search can be examined visually by opening the appropriate .t21 file using ADFMOVIE (SCM:MOVIE).

If it is possible to simply select Transition State Search and set the calculation going, why have we made so much of approximate reaction coordinates to find transition states especially since the direct method finds a real TS while an approximate reaction coordinate finds an approximation to the TS? The reason is that direct TS searches are actually quite difficult to deal with and often fail.

TS are difficult to find. A TS search has a lot in common with a geometry optimization. It is performed stepwise, covergence criteria must be specified and the energy is minimized with respect to the geometry change in most directions. The difference is the energy must be maximized in one direction. One problem is that the TS that you are looking for is not always the one that you find. In some cases nothing is found at all. None of the method for finding TSs available right now is completely satisfactory. Developing new methods of finding TSs is still an active area of research.

The search method that ADF uses often has trouble finding a TS. It works best if the starting geometry is as close as possible to the structure of the TS. This is where the approximate reaction coordinate method comes in. The approximate TS found using an approximate reaction coordinate is often a good starting guess for an direct TS search.

Further Information: For details on some of the algorithms for finding transition states that are commonly used today look at http://theory.cm.utexas.edu/henkelman/research/.
Exercise 5.8: Starting from your best approximation to the TS of the HCN → CNH reaction from the use of an approximate reaction coordinate, run a direct TS search. Compare the geometry and energy that you obtain from the direct method with what you got from the approximate reaction coordinate method. Repeat this calculation except this time take the structure of either the reactants or the products as you starting geometry. Repeat it again with a starting structure that has an H-C-N angle of 120 degrees. Compare the three TS search calculations. If any of these searches runs out of steps before it converges, modify the maximum number allowed through the Detail:Task: Transition State option to allow 100 steps and repeat. If it still fails after 100 steps leave it at that.

5.9 Project: Ethene Polymerization via a Homogeneous Catalyst

The propogation step of ethene polymerization with a homogeneous catalyst is generally believed to have the following mechanism:
Figs/Lab5/Proj_5_2.jpg

where M is usually a transition metal, R is the growing polymer chain, R' and R" are organic ligands and the ethene is bound to the metal through its π-bond.

We will consider the reaction where the metal M is Zr in the +4 oxidation state, the polymer chain is methyl and the R groups are a chlorine atom. Calculate the energy of this reaction and the height of the reaction barrier. To do this you will need to calculate the energy of the reactants ([MeZrCl₂]⁺ with a complexed ethene molecule) and products ([PrZrCl₂]⁺) and the energy and geometry of the transition state. You are free to find the transition state with any method that you like but are advised to define an approximate reaction coordinate and use that coordinate in a linear transit calculation to obtain a good guess of what the transition state looks like. From that good guess a direct transitions state search should be successful.

Use the BP functional and a DZP basis set.

Once you have found the reactants, products and transition state answer the following questions:

What reaction energy do you get? What barrier height?
In the reactants, describe the orientation of the ethene molecule bonded to the Zr atom relative to the methyl group.
In the products, describe the conformation of the polymer chain that is now two carbon atoms longer.
Describe the gometry of the transition state in comparison to the geometry of the reactants and products.
Is your transition state geometry consistent with the proposed mechanism?
Sometimes in reactions of this type the transition states includes a hydrogen atom from the growing polymer chain going very close to the metal atom and forming an "agostic interaction." Do you see this in your transition state? If so, describe what happens.
Another similar catalyst has a barrier of 50 kJ/mol for the ethene insertion reaction. Is this catalyst better or worse than your catalyst?
The polymerization reaction can be terminated in a few different ways. What must be the minimum height of the barriers of these termination reactions so that the polymer produced by this catalyst will have an average molecular weight of at least 100,000 amu? Assume that the polymerization reaction is taking place at 300 K.

File translated from T_EX by T_TH, version 3.85.
On 29 Sep 2009, 15:38.

Lab 5 Chemical Reactions