CP3MD

Molecular Reaction Dynamics: Applications to Triatomic systems

Based on a script revised in Feb 2015 by Pietro Aronica and Dr Paul Wilde. Updated for 2016 by Tristan Mackenzie, Nathan Fitzpatrick, Dr João Malhado and Prof Michael Bearpark

The software

Instructions on how to run the software and obtain the data files required for this experiment are here: Running MD code in MATLAB.

Chemistry computer rooms have been booked for the 2nd year computational exercises when they are running this term, and demonstrators will be available on the days these are timetabled for.

Introduction

The write-up Required

A data sheet template/guide to the write-up and points for discussion required is available here.

Enter data and include screenshots in the data sheet at points indicated.

Submit your completed data sheet via Blackboard by 6 pm on the Friday in the week you carry out the exercise.

Background reading

Note that it is possible to start this exercise and run some trajectories very easily but reference to literature will be needed to understand the concepts introduced.

Atkins and de Paula: Physical Chemistry, 7^th Edition. Topic titles are provided so that you can find the material in later editions of this textbook or in other textbooks.

• Chapter 25: The Rates of Chemical Reactions
• Section 25.5 - T dependence of reaction rates, Arrhenius Parameters their determination and interpretation.
• Chapter 27 Molecular Reaction Dynamics.
• Sections 27.4 Activated Complex Theory, The Eyring Equation.
• 27.5 Thermodynamic aspects of activated complex theory
• 27.6 Dynamics of Molecular Collisions.
• 27.7 Potential Energy Surfaces (PESs)
• 27.8 Some results from experiments and calculations with PESs 

Note that many other physical chemistry textbooks will also cover this material. Editions of Atkins will vary in their coverage.

Objectives

The objectives of this exercise are to study the reactivity of triatomic systems, where an atom and a diatomic molecule collide, through calculating Molecular Dynamics trajectories.

This exercise is designed to make you think about transition states, reaction coordinates and potential energy surfaces. You should discover that, for a chemical reaction to take place, having sufficient energy to overcome a barrier is not enough: the energy must also be in the right vibrational modes at the right time. You'll see trajectories cross the barrier, as if a reaction is going to take place, but then go back to reactants, with energy exchanged between translation and vibration.

The gas phase collision and reaction between an atom and a linear diatomic molecule will be studied. If the reaction proceeds it will form a new diatomic molecule and atom. This exercise builds a level of sophistication upon previous kinetic treatments that you have seen by examining other types of energy (e.g. vibrational) that may be needed for reaction or present in products and also by illustrating how certain conditions lead to reaction profiles which simpler models do not predict.

Relevance to other courses

A lot of the background material for this exercise has been covered in courses you have taken in year 1.

• Molecular Driving Forces (intermolecular interactions - Dr. Ian Gould)
• Chemical Kinetics (collision theory - Dr Oscar Ces)
• Spectroscopy, Quantum Chemistry (vibrational energy, harmonic oscillator - Prof Klug, Dr Wilde).
• Introduction to Physical Organic Chemistry  (Transition states, activation enthalpies, reaction energy profiles - Prof McCulloch).

Having completed the exercise you will be better prepared to understand material in year 3 where the Molecular Reaction Dynamics course (Prof Yaliraki) in particular builds on the ideas and concepts developed here.

Theory

Introduction

It is possible to simulate a chemical reaction by describing the relative motion of the atoms that occurs while the reaction takes place. Since the mass of the atoms is relatively large (compare with the mass of the electron for example), for many chemical reactions it is possible to assume that the motion of the atoms follow Newton's equations of motion, i.e. obey classical mechanics. This is an approximation, since it neglects the fact that molecular vibrations are quantized or quantum tunneling effects, but it is nevertheless a good approximation for most chemical reactions.

The relative motion of the atoms will depend on the interaction between them, which in turn depends on their relative position (for example, how far their nuclei are from each other). Interatomic interactions, as part of a molecule or not, are expressed as a potential energy surface V(r) which represent the potential energy of the system as a function of the atoms' relative positions defined by the general vector r. (For a diatomic molecule V(r) is just the usual potential energy curve that expresses the change in the energy as the two nuclei are displaced from their equilibrium positions by changing the interatomic distance – making the separation too small leads to a rapid rise in potential energy and increasing the bond length also leads to an increase in energy and eventually to dissociation).

The force (variation of momentum p) acting on a given interatomic coordinate r_i will depend on the derivative of the potential energy surface with respect to that coordinate.

${dp_i \over dt} = - { \partial V(r_1,r_2,...)\over \partial r_i}$

By solving the equations of motion we are able to determine the trajectory of the system, i.e. determine the relative position of the atoms at each instant in time r(t). The trajectory represents a path across the potential energy surface (you will see examples shortly), and can be used to gain insight on how the reaction might proceed.

Such simulations have many applications in chemistry and biology and you may well encounter these in higher level lecture courses in year 4.

Note that the program does all these calculations for you. You input bond distances and momenta for the atoms specified and the potential energy surface and the trajectory the system follows across the surface is calculated for you. Your job is to interpret the results.

Atom-Diatom Collisions

In this exercise we will apply this method to study atom diatom collisions. A typical example is depicted below.

Figure 1: An example of a triatomic exchange reaction

In simple terms, the atom m₁ collides with the molecule and forms a new molecule with m₂, while m₃ is detached as a separate atom. 

If we think about how the distances r₁ and r₂ change during this process we can make several deductions. 

Before the collision/reaction:

• r₂ will be approximately constant (it will oscillate slightly if the molecule has vibrational energy).
• r₁ is decreasing steadily as m₁ gets closer to m₂.

After the collision/reaction:

• r₁ will now be approximately constant (but it will oscillate slightly if the molecule has vibrational energy).
• r₂ is increasing steadily as m₃ moves away from the molecule formed from m₂  and m_1.

This means that instead of our simple picture of energy versus reaction coordinate we can map the progress of this reaction on a potential energy surface where potential energy is plotted as a function of both r₁ and r₂. We thus obtain a contour plot like the one shown in Figure 2. For a linear triatomic exchange reaction, the potential function V is a function of only two parameters r₁ and r₂. Thus we can plot the potential energy as a contour plot of V(r₁,r₂) versus r₁ and r₂. A trajectory {r₁(t) r₂(t)} can then be drawn on this potential energy surface.

Figure 2: An example of a reactive trajectory for H + H₂ potential surface

In Figure 2, the small boxes indicate the reactants (r₁ large r₂ small), a transition structure where r₁=r₂, and the products (r₁ small r₂ large). (Note that the variation of potential energy along the combined r₂ and r₁ axes is the typical curve you have seen before). The general triatomic system depicted means we started with atom A and molecule BC and ended with molecule AB and atom C. A reactive trajectory (it’s reactive because it proceeds all the way from reactants to product; not all trajectories will do this) that passes through the transition structure is shown as a wavy line. (The transition structure is a saddle point in the potential energy surface). Notice that the trajectory shown involves vibration in r₂ (the molecule BC has vibrational energy) as the inter-fragment distance r₁ is decreased in the entrance channel (at the top left of the figure) and vibration in r₁ (the product AB has vibrational energy) as the inter-fragment distance r₂ is increased in the exit channel (at the right of the figure). An unreactive trajectory would bounce off the barrier or its surroundings and regenerate the reactants.   

For a given potential surface, the outcome of a dynamics simulation (i.e. the positions r₁(t), r₂(t) and the momenta p₁(t), p₂(t) , p_i =mv_i) at some time t are determined by the initial conditions at time t=0 (i.e. r₁(0), r₂(0) and the momenta p₁(0), p₂(0)).

Only certain initial conditions will lead to a trajectory that passes via the transition structure to the products even, if one has enough energy to surmount the barrier.

In the program you will use, you will run trajectories on different potential surfaces - corresponding to different triatomic chemical systems - for a variety of input conditions.  You will be able to animate the trajectory on the potential surface and examine the distribution of vibrational and translational energy in the products. 

As a demonstration you will first study the potential surface for an atom of H colliding with a molecule of hydrogen to illustrate the principles and to learn how to use the main features of the program. You will then experiment with some more subtle examples.

EXERCISE 1: H + H₂ system

The instructions for this exercise are rather verbose. The objective is to illustrate the various features of the program and then apply them to studying some aspects of this reaction which you report upon.

The first stages of this exercise are designed to familiarise you with the programme. They do not form part of the assessed work directly, but there are questions on the data sheet / writeup template to answer. Please do not rush through this part of the exercise, as it is designed to develop your understanding of potential energy surfaces, reaction coordinates, transition structures, and to introduce the difference between kinetics and dynamics which can be studied experimentally.

Viewing the potential surface

Before you begin experimenting with the different types of collision, note that the script makes reference to r₁ and r₂ whereas in the GUI the atoms are currently labelled A B and C and the distances are r_AB and r_BC. Since you run some trajectories in both forward and reverse directions it is hard to develop a consistent notation. Thus before you run a trajectory, check carefully what you have set up. Your initial molecule will be defined by which of the AB or BC distances you have set to be the bond distance. One distance will be short (the bond distance) and one will be long. 

Run a reactive trajectory

In this part you are given a set of conditions such that the reaction depicted below will proceed. The example takes you through how to use the different options available (the different types of plot) to see what is happening for a particular reaction. Thus you explore the trajectory by looking at contour (2D) and surface (3D) plots and how internuclear distances, potential and kinetic energy and momenta change during the reaction.

The default settings of the positions and momenta are initialised to produce a reactive trajectory (like the one in Figure 2). In other words, these conditions should lead to the reaction below taking place. As you will see in later sections, not all trajectories you run will be reactive.

Figure 3: Initial conditions that lead to reaction

The file you downloaded (lepsgui.m) has a set of initial conditions and will run with those when you first use the command and should show you a contour plot.

To run the example above, first check the initial conditions you have in the GUI. If the initial conditions you have are not as noted above, then change as needed. 

For example in Figure 3 you can see both momenta are -2.7 so if you have r_BC set at 0.74 make sure BC momentum is changed to 0.0.

• Check the initial conditions that are shown in the menu are as follows: (r₁(0) = 0.74  r₂(0) = 2.30  and the momenta p₁(0)=0  p₂(0) = -2.7). 
• A negative value of the momentum corresponds to a velocity that decreases the interatomic distance.  Thus in this case we are starting the trajectory by giving the system some velocity in the direction that decreases r₂ (refer to the diagram above).

You can see from a comparison of this data with Figure 3 that in this case BC is the molecule with a bond distance = 0.74 and A is the atom that collides with it. Thus here r₁ = r_BC.

• The program updates the image on the right whenever you press the “Update” button at the top (see instructions). When you started the program, it will have automatically calculated the trajectory with the initial conditions outlined above and shown you the contour plot. To see other plots, such as plots of velocities and positions against time, or a 3D surface plot of the reaction, select that option from the drop down menu at the bottom left and press update. You can have a better look at the surface plot by rotating it (the button on the bar at the top with the arrow circling the cube). Spend some time to have a look at the various options.

• Select “Surface Plot”, which gives you a plot of the Potential Energy Surface (PES). The thick black line in the middle is the path of the reaction. Try to identify the reactant and product channels. Observe how the line is wavy: this indicates that the diatomic molecule is vibrating. The most important point to observe is conservation of energy: the system will always oscillate between energy contours with the same value.

• To animate the geometry, select the last option, “Animation”, in the drop down menu at the bottom left. You will see that initially r₂ decreases as H₃ approaches H₁-H₂, then after the transition state, H₁ leaves and H₂-H₃ is formed in a vibrationally excited state.

• This same information can be displayed in a different form by selecting “Internuclear Distances vs Time” in the same menu. In this graphical display the values of r₁(t), r₂(t) are given (Y axis) against time t (X axis). Notice that r₂ (r_AB here) decreases from the original value of 2.3 during the first 0.38 time units, then it has an oscillatory behaviour for the remainder of the trajectory corresponding to H₂-H₃ vibration.  In contrast r₁ (r_BC here) has a slight oscillatory behaviour initially, then after the transition state, it grows in value as H₁ leaves. Notice that r₁(t) = r₂(t) at ≈ 0.38 s. This is effectively the transition state.

Figure 4: r₁(t) and r₂(t) are given (Y axis ) against time t (X axis). Note that only two distances are shown on this plot. The GUI shows three.

• Next, select “Potential Energy vs Time” in the same menu.  It will show an oscillatory behaviour rising in energy until the transition state is reached. The potential energy then falls as the trajectory moves into the product / exit channel. The energy is continually switching from potential energy to kinetic energy during the trajectory. At the turning point of a trajectory one is at the maxima of the oscillations. At the bottom of the trough of the oscillations the potential energy is low and the molecule is moving faster so the energy has gone into kinetic energy. You can see this with “Kinetic Energy vs Time”. In this case the kinetic energy reaches its minimum value at the transition state (i.e. the system is moving slowly).
• Finally the momenta p₁(t) and p₂(t) (p_i =mv_i),  can be displayed using “Internuclear Momenta vs Time”. This result is not easy to interpret because the momenta correspond to bond stretching displacements r₁ and r₂ as illustrated below.

Figure 5: Momenta p₁ and p₂ correspond to those that change the interatomic distances, not the atoms themselves

Thus p₁ is the momentum in the coordinate r₁. Notice that the momentum in the atom 2 results from both momenta p₁ and p₂.   Thus if we want to collide atom 3 with 1-2 (by giving some p₂) in such a way that the diatom 1-2 does not vibrate, we must give some momentum in p₁ as well.  Following the momentum distribution from the transition state provides some insight as we now discuss.

The important observations are as follows:

Initially p₁(0)=0 and p₂(0)=-2.7.  At the beginning of the trajectory p₂(t) decreases and p₁(t) takes on an oscillatory behaviour. p₁(t) oscillates because the initial momentum in p₂ also forces r₁ to change as shown by the picture above. After the transition state at t=0.38, p₁(t) rises rapidly and remains constant (corresponding to translation) as r₁ becomes large as H₁ leaves.  In contrast p₂ initially corresponds to translation. After the transition state p₂ oscillates (corresponding to vibration).           

Dynamics from the transition state region

Figure 6: Gradients at the TS region

The transition state is defined as the maximum on the minimum energy path linking reactants and the products. This point on the potential energy surface has the special property that dV(r)/dr (the gradient of the potential) is zero, and the energy goes down most steeply along the minimum energy path linking reactants and products. Consequently, if one starts a trajectory exactly at the transition state, with no initial momentum, it will remain there forever.  However, if one changes the geometry by a small amount in the direction of the products it will roll towards the products (and similarly for the reactants).  One way of locating the transition state is to start trajectories near the transition state and see whether they "roll" towards the reactants or products. The reaction path (minimum energy path or mep) is a very special trajectory that corresponds to infinitely slow motion (i.e. the velocity always reset to zero). Once the transition state has been located, one may run the very special trajectory that corresponds to the mep.

Since the H + H₂ surface is symmetric, the transition state must have r₁  = r₂.  If we start a trajectory on the ridge r₁ = r₂ there is no gradient in the direction at right angles to the ridge, thus the trajectory will oscillate on the ridge and never fall off.  As we will see, this fact can be used to locate the TS geometry. If we start a trajectory at r₁ = r_ts+δ and r₂ = r_ts the trajectory “falls off” the r₁ = r₂ ridge. We begin by running trajectories from near the TS geometry.

Discussion: What value does the total gradient of the potential energy surface have at a minimum and at a transition structure? Briefly explain how minima and transition structures can be distinguished using the curvature of the potential energy surface.

Up to this point, you have been learning how to use the program to visualise a triatomic reaction. From now on you will be exploring the factors that influence whether or not a trajectory is reactive and why. You should make a note of the key data and results because they should be included in your report. From this point on the work you do should be included in your report and will be assessed.

Trajectories from r₁ = r₂: locating the transition state

• Use the menus and buttons on the left of the GUI to modify the initial conditions.
• First set the momenta to exactly 0. Type in values of the momenta corresponding to p₁, p₂. You should supply exactly 0.0 here.
• Then type in values of r₁, r₂.  You should supply one choice of values of r1 and r₂ that are exactly equal in the region of 0.85 to 0.95. To see the reaction path, select “Contour” or “Surface Plot” on the bottom menu, then press “Update” to visualise. If you select “Animation”, you will see that the system undergoes a periodic symmetric vibration of the form shown in Figure 7.

Figure 7: Symmetric vibration about the TS

• Now select “Internuclear Distances vs Time”. You will see that the plots for r₁ and r₂ are superimposed because they are oscillating exactly in phase. You will thus only be able to see two of the three colours on this example. The geometry about which they are oscillating is the transition state geometry. Determine this value of r₁ = r₂ as accurately as you can from the graph. As further proof that you have found the value, try typing it in the initial conditions: it should give you a straight line without oscillation. We shall refer to this value as r_ts. Make a note of this value on the data sheet. 

Trajectories from r₁ =  r_ts+δ, r₂ = r_ts

• Now we will run trajectories in the forward and reverse directions from r₁ = r₂ = r_ts with p₁ = p₂ = 0.  This can be accomplished by running trajectories with r₁ = r_ts and r₂ = r_ts+0.01. This should terminate at r₂ > 2.0 with r₁ oscillating i.e. at H₁-H₂ + H₃. Convince yourself that a trajectory with r₁ = r_ts+0.01 and r₂ = r_ts will give similar results terminating at r₁ > 2.0 with r₂ oscillating i.e. at H₁+ H₂-H₃.

Use different plots (Contour, Surface, Animation) to analyse the system.

• Finally, look at “Internuclear Distances vs Time” and note the final values of the positions r₁(t) r₂(t) and the average momenta  p₁(t) p₂(t) (using “Internuclear Momenta vs Time” at large t for later use. Write these down on the data sheet.

 Consider the results of the trajectory with r₁ = r_ts and r2 = r_ts+0.01 which terminated at r₂ > 2.0 with r₁ oscillating i.e. at H₁-H₂ + H₃ Thus asymptotically, the r₁ oscillates about the H₁-H₂ equilibrium bond length (0.74) as r₂ became large. The momentum p₁ oscillates around 1.25 (i.e. between 0.8 and 1.5) while p₂ becomes constant at around 2.5. By simply reversing the signs of the final momenta we can run the trajectory backwards to give a reactive trajectory. This data tells us that any trajectory will be reactive for r₁ = 0.74, r₂ = 2.0, with -1.0 < p₁ < -1.5 and p₂ = -2.5

The mep from r₁ = r₂

Next, we will look at the mep. Leave the geometry and momenta as in the previous exercise (i.e. the positions r₁ = r_ts+.01,  r2 = r_ts and the momenta p₁=0 and p₂=0), then change the calculation type from dynamics to MEP. You will see that the trajectory simply follows the valley floor to H₁+ H₂-H₃

Generating reactive trajectories 

How do we determine the conditions for a reactive trajectory that starts in the region of the reactants and passes near the TS region?  As you will discover, simply adjusting the collision parameters so that the system has enough energy to reach the TS is not enough. The energy must be distributed correctly. 

• Change the computation type back to Dynamics.

• From our computations that began at the transition state with r₁ = r_ts and r₂ = r_ts+ 0.01 we know that a reactive trajectory will correspond to r₁ = 0.74, r₂ = 2.0, with -0.8 < p₁ < -1.5 and  p₂ = -2.5
• Run several trajectories with fixed values of r₁ = 0.74, r₂ = 2.0, p₂ = -2.5 but with -0.8 < p₁ < -1.5.  They should all be reactive. Note that p₁ is varied here. Make sure your conditions are correct before you start. Report on the data sheet conditions that you used and the characteristics of the trajectories that you observe.

• Run the trajectory with r₁ = 0.74   r₂ = 2.0, p₁ = -1.25 and p₂ = -2.5. What do you note about the behaviour in the entrance channel? Report on the data sheet. 

• Run a trajectory with r₁ = 0.74, r₂ = 2.0, p₁ = -1.5, p₂ = -2.0. Report on the data sheet. Is the trajectory reactive or unreactive? Why? 

• Run a trajectory with r₁ = 0.74, r₂ = 2.0, p₁ = 1.5, p₂ = -2.5. Report on the data sheet. Is the trajectory reactive or unreactive? Why?  

• Run a trajectory with r₁ = 0.74, r₂ = 2.0, p₁= -2.5, p₂ = -5.0 (i.e. similar to the trajectory where we force pure translation in the input channel but with the momenta twice as large). The result is remarkable:  The trajectory actually crosses the transition state region but does not go on to products even though we have given it more than twice the momentum needed to surmount the barrier. It is instructive to use “Animation” for this example.  You should see the bond in the product actually forms but then the system reverts back to the reactants. Report on the data sheet and describe what happens. 

Figure 8: Trajectory that re-crosses the TS region

• Running a trajectory with r₁ = 0.74, r₂ = 2.0, p₁ = -2.5, p₂ = -5.2 produces an unusual result. Report on the data sheet. 

Discussion: In the transition state theory of reaction rates it is assumed once the system reaches the transition structure it goes on to produce products.  However transition state theory usually overestimates the reaction rate. Why?

Guidance on your report for Part 1: the H and H₂ system

Discussion of the first part of the exercise

The reaction between H and H₂ has been studied here. Use the data sheet and the observations that you recorded there to structure your results and discussion. Explain why a particular outcome is observed for a particular set of conditions.

Your focus should be to report the outcomes of the trajectories that you have run, whether they are reactive or unreactive and any key observations that particular trajectories show.  

EXERCISE 2: F + H₂ system

• Select the F + H₂ potential surface by changing the atom types in the GUI. View the surface plot (potential surface). You need not worry about the conditions as the point is to view the surface, not a trajectory.
• You will see that the reaction is highly exothermic (i.e. the exit channel is much lower in energy than the entrance channel since the HF bond is much stronger than the HH bond). Furthermore, there appears to be no transition state.
• Because the entrance and exit channels are so different, you may need to experiment with the 3D view rotation (see instructions for part 1) to find the most informative point of view.

The mep has the following form:

Close inspection shows that there is a barrier in the entrance channel (and consequently a transition state) but the barrier height is so small that it is not within the resolution of the program you are using. At the transition state r₁ (F-H₂) = 1.80 and r₂ (FH-H) = 0.75. Because the reaction is exothermic, the transition state is reactant-like and occurs in the entrance channel (an early barrier).

With an exothermic reaction such as this one you may find it useful to explore the potential energy surface and the trajectories using the data cursor. Figure 9 below shows the icon on the toolbar. Clicking on the icon allows you to interrogate part of the surface to find coordinates and energy values (indicated by level – see the Figure).

Figure 9: The data cursor allows you to map the surface and find coordinates (r_AB, r_BC and energy (“level”).)

Run an mep computation from the TS

Set the momenta to p₁ = 0 and p₂ = 0 and offset the distances from the TS to r₁ = 1.78, r₂ = 0.75. Then select the calculation type MEP and run the trajectory in the usual way. You will see that the trajectory simply follows the valley floor to  F-H + H. Note: You will need a large number of steps (>10000) to reach the valley floor.

Include a brief discussion and illustration of the mep in your report.

Dynamics computation from the TS: Release of the energy of exothermicity

• Change the calculation type back to Dynamics.
• Run the trajectory again. Examine the results with “Animation”, and look at the “Internuclear Distances vs Time” and “Internuclear Momenta vs Time”.
• What do you conclude about the release of energy in this reaction?

Include a brief discussion of the energy release for this exothermic reaction with suitable illustrations in your report.