Molecular Reaction Dynamics

In this lab we are using molecular dynamic trajectories to study the reactivity of triatomic systems, namely when an atom collides with a diatomic molecule.

ADD CONTENTS TABLE

EXERCISE 1: H + H₂system

For this exercise, the following reaction is analysed and discussed: H_A + H_BH_C→ H_AH_B + H_C, where r₁ = BC and r₂ = AB.

Dynamics from the transition state region

Question 1: On a potential energy surface diagram, how is the transition state mathematically defined? How can the transition state be identified, and how can it be distinguished from a local minimum of the potential energy surface?

Answer 1: On a potential energy surface diagram the transition state is defined, mathematically, as a saddle point. A transition state is the maximum energy present in the minimum energy pathway taken from reactants to products. A saddle point is of symmetric configuration at which the distance between three separate atoms is the same.

A saddle point is different to a local minima. For both, the gradient of the potential energy surface is zero, is given as the first partial derivative:

∂V(r_i)/ ∂r_i = 0

However, the saddle point and local minima may be distinguished though inspection of their second partial derivatives. For a local minima, calculating partial second derivatives in orthogonal directions to the saddle point will give a result that is greater than 0, independent from which plane it has been measured. By contrast, for the saddle point, is second partial derivative will either be greater than zero or smaller than zero dependent on which plane it has been measured. This explanation is summarised below:

For a local minima: ∂²V(r_i)/ ∂r_i² > 0

For a saddle point: ∂²V(r_i)/ ∂r_i² > 0 and ∂²V(r_i)/ ∂r_i² < 0

An example of a saddle point can be seen in the figure below:

Figure 1: An example of a saddle point, orthogonal directions from which have opposite signs for their partial second derivative.

Question 2: Report your best estimate of the transition state position (r_ts) and explain your reasoning illustrating it with a “Internuclear Distances vs Time” plot for a relevant trajectory.

Answer 2: The transition state occurs at the saddle point, which is of symmetric configuration. Hence, r₁ = r₂as the H + H₂ surface is symmetric, where, r₁= BC and r₂ = AB. At the transition state p₁ = p₂ = 0 g.mol^-1.pm.fs^-1. Through trial and error an estimate for the transition state position is when AB = BC = 90.75 pm. Using this distance produces the contour plot below. It can be seen that there is no gradient at right angles to the point of the transition state.

Figure 2: A contour plot showing the transition state. There are no oscillations or directions of motion towards either reactants or products, thus the system is stationary

The plot of 'Internuclear Distances vs Time' below shows the distances of A-B, B-C, and A-C constant with time, consolidating the fact that this position is the transition state position. Running an animation with the r₁ = r₂= 90.75 pm shows the species as stationary, as there are no fluctuations of interatomic distance with time.

Figure 3: A plot of Internuclear Distances vs Time

Question 3: Comment on how the mep and the trajectory you just calculated differ.

Answer 3: The initial conditions are set so that the system is slightly displaced from the transition state and so that there is zero initial momenta, where r₁ = r_ts + 1 pm, r₂ = r_ts and the momenta p₁ = p₂ = 0 g.mol^-1.pm.fs^-1. This can be shown in the figures below, where r₁ (BC) = 91.75 pm and r₂ (AB) = 90.75 pm. To obtain a complete MEP the step number was increased to 3000.

Figure 4a: Contour Plot (Calculation Type: Dynamics)		Figure 5a: Contour Plot (Calculation Type: MEP)

Figure 4b: Internuclear Distance vs Time Plot (Calculation Type: Dynamics)		Figure 5b: Internuclear Distance vs Time Plot (Calculation Type: MEP)

Figure 4c: Momenta vs Time Plot (Calculation Type: Dynamics)		Figure 5c: Momenta vs Time Plot (Calculation Type: MEP)

It can be deduced that the trajectory is greatly dependent on calculation type. Using the MEP calculation, the momenta are reset to zero after each step. This provides an unrealistic account of the motion of the atoms during the reaction. This explains the absence of any type of oscillation in the MEP calculation, Figure 5a, compared to the Dynamics calculation figure. The Dynamics calculation accounts for the vibrational energies possessed by the molecules, hence explaining the presence of oscillations along the trajectory in Figure 4a.

Plots of "Internuclear Distances vs Time” and “Momenta vs Time” are also provided and more clearly depict the differences in the calculation types.

The Dynamics calculation accounts for the extra kinetic energy obtained after bond formation, which manifests itself into excess vibrational energy. This can be seen in the small oscillations of the A-B bond in figure 4b. Using the MEP calculation, this A-B bond distance remain constant over time, indicating this excess vibrational energy has not been accounted for. The increased momentum of the B-C bond is indicative of H_C moving with translational energy and the oscillating momentum of the A-B bond is indicative of the new H_AH_B bond moving with vibrational energy.

Figure 5c shows a straight line at 0. The momenta is zero for the reaction pathway. This calculation type resets the momenta to zero after every step hence this result is expected. Figure 4c shows an oscillation of momentum with time, as through this calculation type, momenta is not reset to zero after each step.

If the initial conditions were changed so that r₁ = r_ts pm and r₂ = r_ts + 1pm, the blue and orange lines would swap positions in figures 4b, 4c, 5b, and 5c.

Reactive and unreactive trajectories

Question 4: Complete the table below by adding the total energy, whether the trajectory is reactive or unreactive, and provide a plot of the trajectory and a small description for what happens along the trajectory. What can you conclude from the table?

With the initial positions r₁ = 74 pm and r₂ = 200 pm, trajectories were run with the following momenta combination:

p₁/ g.mol^-1.pm.fs^-1	p₂/ g.mol^-1.pm.fs^-1	E_tot/ kJ.mol^-1	Reactive?	Description of the dynamics
-2.56	-5.1	-414.280	Yes	The kinetic energy of the system is sufficient to over come the activation barrier and hence the reaction proceeds. Originally, the reactants have little vibrational energy, reflected by the lack of oscillatory behaviour. After the reaction occurs, the products have more oscillatory behaviour, meaning that they have greater vibrational energy. The BC bond is broken and the AB bond is formed.
-3.1	-4.1	-420.077	No	H_A approaches an oscillating H_BH_C molecule, however the activation energy exceeds the kinetic energy and thus the reaction does not occur. Both reactants remain unchanged.
-3.1	-5.1	-413.977	Yes	This trajectory is similar to that of row 1, with the only difference being that the reactants possess more vibrational energy, hence the presence of oscillatory behaviour as depicted in the figure.
-5.1	-10.1	-357.277	No
-5.1	-10.6	-349.477	Yes