File:Q9 2 yc16318.png

2020-05-08T20:36:08Z

Yc16318:

2020-05-08T19:33:56Z

Yc16318: /* Molecular Reaction Dynamics */

=Molecular Reaction Dynamics=
In this report, the the molecular reaction dynamics for H+H2 and F-H-H systems will be investigated in detail. The physical properties of the systems, such as its transition state and activation energy, will be estimated through running dynamic and MEP calculations.

== Exercise I ==
===Q1===
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?

The transition state is mathematically defined as the saddle point, which is a point in a multivariate plane where the second derivative along all coordinates, except for one, is greater than zero. On the potential surface, the reaction trajectories go through a pass of local minima, which includes the saddle point. Therefore, it can also be seen as the maximum minimum energy of the local minima.

[[File:Saddle_yc16318.jpg ‎|thumb|center|A saddle point [[https://www.math.tamu.edu/~tkiffe/calc3/saddle/saddle.html]]]]

For all critical points such as the local maxima, minima and saddle points, all the first partial derivatives will be equal to zero (i.e. the gradient of the point is zero). All second derivatives for local minima will be positive, all second derivatives for local maxima will be negative. However, for the transition state, some of its second derivatives will be positive, while the rest will be negative.

The local minima and the saddle point can be determined using the Hessian matrix:

[[File:Saddle_yc16318.jpg ‎|thumb|center|Hessian Matrix [[https://math.stackexchange.com/questions/2415123/maxima-minima-and-saddle-points]]]]

Hessian_matrix_yc16318.png

If we let D1 = fxx and D2 = det H, we can distinguish between the local minima, local maxima, and saddle point[[https://math.stackexchange.com/questions/2415123/maxima-minima-and-saddle-points]]:
1) If D1(a,b) >0 and D2(a,b) >0, then f has a local minimum at (a,b).
2) If D1(a,b) <0 and D2(a,b) >0, then f has a local maxima at (a,b).
3) If D2(a,b) <0, then f has a saddle point at (a,b)

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

rts≈ 90.8 pm is the best estimate.

The transition state can be estimated by looking at the Internuclear Distances vs Time for when the trajectory oscillates on the ridge only (i.e. r1 = r2) 

[[File:TS1_yc16318.png]]

We can see that A-B overlaps with B-C due to the system's symmetric vibration, and that A-C=2 A-B. The internuclear distances are oscillating. This is because that as the atoms vibrate, they stretch and contract from the transition state position. Therefore for the second estimation, the transition state position rts can be found by finding the average value of the peaks and the troughs of A-B and B-C, which is (152 + 70)/2 = 111. However, 111 is still far from the transition state because we can still see relatively strong oscillations, suggesting it is not at the saddle point where it can remain forever.

[[File:Surface_Plot_TS2.png]]

This is because that our estimation is inaccurate as the potential curve follows the shape of an asymmetric oscillator for large displacements. The lowest point of the curve cannot be found be simple average at this point.

[[File:Surface_Plot_TS.png]]

Using the contour plot, we can estimate r1 using a point inside the contour line (92 pm). In this way, the oscillation can be approximated using harmonic oscillator, where the lowest point can be estimated by the average of the peaks of the troughs.

[[File:IT_TS.png]]

Therefore, r1 can be estimated to be (92 + 89.6)/2=90.8, which gives an Internuclear Distances vs Time plot with straight lines.

[[File:Animation_TS_90_8.png]]

This suggests that one is no longer oscillating and sits on the saddle point. The transition state is estimated.

===Q3===
 Comment on how the mep and the trajectory you just calculated differ. 

[[File:Mep_1.png|thumb|center|Mep calculation]]

In the MEP calculation, inertial motions are disregarded. As we can see in the diagram, B-C decreases to the bonding distance and then remains constant, while A-B increases throughout. There is no oscillation occurring. Physically, the H1 is moving away from H2-H3, and H2-H3 stays constant after it has decreased to around 75 pm. All energy is transnational. 

[[File:Dyn_1_yc16318.png|thumb|center|Dynamic calculation]]

In the dynamic calculation, however, the inertial motions are taken into consideration, which includes the vibration energy into the system. We can see that, as H1 moves away, H2-H3 first decreases, then starts oscillating. Both transnational and vibration energy are included in the system.

Even though the MEP calculation depicts the trajectory with lowest possible energy, it is not realistic. It does not factor in the classical motions of the atoms. However, it gives us further information about the reaction and allows us to estimate how much vibration energy there actually is as simulated in the dynamic calculation. We can infer that, as the system moves slower, it will become more and more similar to the trajectory shown in MEP. Therefore, the slower it moves, the less the hydrogen molecule vibrates.

To gain more in-depth knowledge of the difference between the MEP and the dynamic calculations, we should look at the momentum vs. time plot as well as the internuclear distance vs. time plot for both types of calculations.

[[File:MT_MEP.png|thumb|left|Momentum vs. time, MEP calculation]]

[[File:IDT_MEP1.png|thumb|center|Internuclear distance vs. time, MEP calculation]] 

For the MEP calculations, the momentum is always zero as the motion is assumed to be infinitely small. The rate of increase of AB distance lowers over time and the BC distance first decreases and then stays constant at r2=75 pm. The AB distance approaches 500 pm for large t values. This is because that the potential surface flattens for large AB distance at constant BC distance, hence the velocity of one 'rolling' down the potential surface also decreases. Eventually, the AB distance will also be constant with time.

[[File:MT_DY.png|thumb|left|Momentum vs. time, dynamic calculation]]

[[File:IDT_DY2.png|thumb|center|Internuclear distance vs. time, dynamic calculation]] 

For the dynamic calculation, AB has an increased then constant momentum and BC has an oscillating momentum. This suggest that AB (ie. the speed of H1 moves away from H2-H3) increases at a constant speed once it reaches the 'flattened' potential surface and BC oscillates constantly. Thus, the AB distance becomes linearly proportional to time for large t values.

====Extra====
When the initial conditions and the final conditions are reversed, AB distance and BC distance are swapped, it is expected that the process will be exactly the reverse of the previous scenarios with the plots being the perfect mirror images of each other. 

[[File:Reversed_DT_yc16318.png|thumb|left|Internuclear Distance vs. time, Dynamic calculation, reversed]]

[[File:Reversed_MT_yc16318.png|thumb|center|Momenta vs. time, dynamic calculation, reversed]]

We can see from the plots that, while it resembles the reversed plot of the previous condition, it is not exactly the same. This could be due to the fact that the positions and the momenta of the final conditions were not perfectly estimated. This can also be supported by looking at its contour plot, where the trajectory doesn't pass the transition state due not enough energy in the system.

[[File:Reversed_SP_yc16318.png|thumb|center|Contour plot, dynamic calculation, reversed]]
 
activation energy calculation

===Q4===
Complete the table above 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? 
====Momenta vs. reactivity table====
{| class="wikitable" border="1"
|-
! p1/ g.mol-1.pm.fs-1 !! p2/ g.mol-1.pm.fs-1 !! Etot/ Kj.mol-1 !! Reactive? !! Description of the dynamics !! Illustration of the trajectory 
|-
| -2.56|| -5.1 || -414.280 || Yes || H3 approaches the non-vibrating H1-H2, forming a vibrating H2-H3 as H1 moves away. Products are formed. ||[[File:SP_C1-414_280.png|300px]]
|-
| -3.1 ||-4.1 || -420.077 || No || H3 collides with the noon-vibrating H1-H2. However, it does not move past the transition state as the potential energy exceeds the energy of the system. No products are formed, H3 moves away again. ||[[File:SP_C2-420_077.png|300px]]
|-
| -3.1 ||-5.1 || -413.977 || Yes || Similar to the first condition. H3 approaches the non-vibrating H1-H2, forming a vibrating H2-H3. The vibration is stronger than that in the first condition as the system possesses higher initial momentum. ||[[File:SP_C3-413_977.png|300px]]
|-
| -5.1 ||-10.1 || -357.277 || No || H3 approaches H1-H2. A H2-H3 is formed and vibrates strongly as H1 moves away. However, the system reverts back to the reactants and the H1-H2 is formed again. Eventually, H3 moves away, leaving a vibrating H1-H2 bond. The momentum of the system is large enough to allow the trajectory to cross the transition state twice, but no reaction happens in the end. The vibration is significantly large comparing to the previous conditions due to higher momentum. ||[[File:SP_C4-357_277.png|300px]]
|-
| -5.1 ||-10.6 || -349.477 || Yes ||H3 approaches H1-H2. At first, H2-H3 bond is formed and vibrates strongly as H1 moves away. However, the H1-H2 is not yet fully broken and vibrates once again. Eventually, a vibrating H2-H3 is formed as H1 moves away. The trajectory crosses the transition state three times before the products are formed. ||[[File:SP_C5-349_477.png|300px]]
|}

 

====Conclusion of the table====
From the table, we can conclude that Etot and the momenta of AB and BC are not the only factors for whether a reaction will happen. From the conditions in the table, we can see that, when p1 is kept constant, the condition with bigger p1/p2 difference is reactive, even though it has lower Etot. This could suggest that the difference between p1 and p2 might also play an important role in whether a system is reactive.

===Q5===
Given the results you have obtained, how will Transition State Theory predictions for reaction rate values compare with experimental values? 

In the transition state theory (TST), only the states of the reactants are considered for the rate of the reaction. It is assumed that, once the transition state is passed, the reaction will proceed to exit from the product channel irreversibly. Therefore, only one crossing of the transition state is allowed. Thus, the rate of the reaction can be approximated by the rate of molecules passing through the transition state.

However, as we can see from the previous simulations, it is possible for the reaction trajectories to pass the transition state twice or more if the particles possess large enough momenta. Therefore, we need a sufficiently large t to determine whether a certain pathway is reactive.

Additionally, because we are treating the nuclei as classical point masses, we expect the motions of the particle to be classical even though energy quantisation can be taken into account by the partition functions. Therefore, TST is unable to reflect the quantum effects in a system, such as the tunnelling effect.

London-Eyring-Polanyi-Sato (LEPS) potential 

== Exercise II: F - H - H system ==
===Q6===
By inspecting the potential energy surfaces, classify the F + H2 and H + HF reactions according to their energetics (endothermic or exothermic). How does this relate to the bond strength of the chemical species involved? 

Looking at the surface plot for a FHH system where FH is AB and HH is BC, we can clearly see that the potential surface for when AB distance is small (F-H is bonded) is much lower than that of when BC is small (H-H is bonded). This is because of the big difference in electronegativity for the F and H atoms, which are 3.98 and 2.20, respectively[[https://www.webelements.com/]]. Therefore, the F-H bond is strongly ionic in nature. Such as strong bond possesses much higher bond dissociation energy compare to the covalent H-H bond (H-F = 567 kJ mol-1 and H-H = 436 Kj mol-1)[[https://pubs.acs.org/doi/abs/10.1021/jp060149i]], leading to a much lower potential energy surface.

[[File:SP_FHH_yc16318.png|thumb|center|A surface plot of the FHH system, AB is FH; BC is HH]] 

Thus, the reaction of F + H2 forms product FH +H, which exits from the side with small AB distance and much lower potential surface, is energetically favourable and exothermic; reversely, the reaction of FH +H to F + H2 requires energy input and is therefore endothermic.

===Q7===
Locate the approximate position of the transition state.

The position of the transition state is: rHF =181.104 pm, rHH =74.501 pm. 

The transition state can be estimated using the contour plot. First the position is estimated to be: rHF = 147 pm, rHH = 76 pm by looking at the contour plot. Then the values can be adjusted by generating contour plot using MEP calculation with 0 momentum for both p1 and p2 as shown below.

[[File:FindingTS_YC16318_1.png|thumb|left|finding TS: in this case, the AB value needs to be increased.]]
[[File:TS_FH_H_yc16318.png|thumb|center|In this case, the AB value needs to be decreased.]]

 
 
In this end, the trajectory becomes static at rHF =181.104 pm, rHH, suggesting that we have found the transition state.
 
 

[[File:TS_DT_YC16318.png|thumb|left|The Internuclear Distance vs. time plot for the transition state]]
[[File:TS_FHH_stationary_yc16318.png|thumb|center|The surface plot of the transition state]]

 
 
===Q8===
Report the activation energy for both reactions.
 

{| class="wikitable" border="1"
|+
! Reaction !! Activation energy / kJ mol-1 !! Transition State Energy / kJ mol-1
|-
| H2 + F || 2.486 || -431.94
|-
| FH + H || 127.055 || -431.168
|}

Because we have found that the transition is where rHF (r1) =181.104 pm, rHH (r2) =74.501 pm, we can use the MEP calculation to find the energy vs. time plot at (r1 = rrt + 1 pm, r2 = rrt) and (r1 = rrt, r2 = rrt + 1 pm). The activation energy of H2+F and FH+H can be found from the δE, the difference in energy between the final and initial state, from the two aforementioned plots respectively.

[[File:MEP_ET_F_HH_yc16318.png|thumb|center|Ea= 2.486 kJ mol-1 can be calculated by subtracting the initial value -431.94 by the final value -434.426]]

[[File:MEP_ET_FH_H_yc16318.png|thumb|center|Ea= 127.055 kJ mol-1 can be calculated by subtracting the initial value -431.168 by the final value -558.223]]

===Q9===
In light of the fact that energy is conserved, discuss the mechanism of release of the reaction energy. Explain how this could be confirmed experimentally.
answer

===Q10===
Discuss how the distribution of energy between different modes (translation and vibration) affect the efficiency of the reaction, and how this is influenced by the position of the transition state.
answer

File:Hessian matrix yc16318.png

2020-05-08T19:33:15Z

Yc16318:

MRD:yc16318

2020-05-08T18:24:39Z

File:SP FHH yc16318.png

2020-05-08T17:52:17Z

Yc16318:

MRD:yc16318

2020-05-08T16:12:26Z

Yc16318: /* Extra */

=Molecular Reaction Dynamics=
In this report, a detailed explanation of the molecular dynamics will be given.

== Exercise I ==
===Q1===
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?

The transition state is mathematically defined as the saddle point. (evidence&graph)

All partial derivatives will be equal to zero for both saddle point and local minima. (i.e. the gradient of the point is zero)

It can be distinguished from a local minimum by by calculation of its second derivative. This is because that, at the saddle point, the second derivative along all the coordinates, except for one, is greater than zero. (details and google, consult yr 1 MPC), hessian matrix and the determinant

can use discriminant to distinguish local minima and saddle point.
https://math.libretexts.org/Bookshelves/Calculus/Map%3A_University_Calculus_(Hass_et_al.)/13%3A_Partial_Derivatives/13.7%3A_Extreme_Values_and_Saddle_Points

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

rts≈ 90.8 pm is the best estimate.

The transition state can be estimated by looking at the Internuclear Distances vs Time for when the trajectory oscillates on the ridge only (i.e. r1 = r2) 

[[File:TS1_yc16318.png]]

We can see that A-B overlaps with B-C due to the system's symmetric vibration, and that A-C=2 A-B. The internuclear distances are oscillating. This is because that as the atoms vibrate, they stretch and contract from the transition state position. Therefore for the second estimation, the transition state position rts can be found by finding the average value of the peaks and the troughs of A-B and B-C, which is (152 + 70)/2 = 111. However, 111 is still far from the transition state because we can still see relatively strong oscillations, suggesting it is not at the saddle point where it can remain forever.

[[File:Surface_Plot_TS2.png]]

This is because that our estimation is inaccurate as the potential curve follows the shape of an asymmetric oscillator for large displacements. The lowest point of the curve cannot be found be simple average at this point.

[[File:Surface_Plot_TS.png]]

Using the contour plot, we can estimate r1 using a point inside the contour line (92 pm). In this way, the oscillation can be approximated using harmonic oscillator, where the lowest point can be estimated by the average of the peaks of the troughs.

[[File:IT_TS.png]]

Therefore, r1 can be estimated to be (92 + 89.6)/2=90.8, which gives an Internuclear Distances vs Time plot with straight lines.

[[File:Animation_TS_90_8.png]]

This suggests that one is no longer oscillating and sits on the saddle point. The transition state is estimated.

===Q3===
 Comment on how the mep and the trajectory you just calculated differ. 

[[File:Mep_1.png|thumb|center|Mep calculation]]

In the MEP calculation, inertial motions are disregarded. As we can see in the diagram, B-C decreases to the bonding distance and then remains constant, while A-B increases throughout. There is no oscillation occurring. Physically, the H1 is moving away from H2-H3, and H2-H3 stays constant after it has decreased to around 75 pm. All energy is transnational. 

[[File:Dyn_1_yc16318.png|thumb|center|Dynamic calculation]]

In the dynamic calculation, however, the inertial motions are taken into consideration, which includes the vibration energy into the system. We can see that, as H1 moves away, H2-H3 first decreases, then starts oscillating. Both transnational and vibration energy are included in the system.

Even though the MEP calculation depicts the trajectory with lowest possible energy, it is not realistic. It does not factor in the classical motions of the atoms. However, it gives us further information about the reaction and allows us to estimate how much vibration energy there actually is as simulated in the dynamic calculation. We can infer that, as the system moves slower, it will become more and more similar to the trajectory shown in MEP. Therefore, the slower it moves, the less the hydrogen molecule vibrates.

To gain more in-depth knowledge of the difference between the MEP and the dynamic calculations, we should look at the momentum vs. time plot as well as the internuclear distance vs. time plot for both types of calculations.

[[File:MT_MEP.png|thumb|left|Momentum vs. time, MEP calculation]]

[[File:IDT_MEP1.png|thumb|center|Internuclear distance vs. time, MEP calculation]] 

For the MEP calculations, the momentum is always zero as the motion is assumed to be infinitely small. The rate of increase of AB distance lowers over time and the BC distance first decreases and then stays constant at r2=75 pm. The AB distance approaches 500 pm for large t values. This is because that the potential surface flattens for large AB distance at constant BC distance, hence the velocity of one 'rolling' down the potential surface also decreases. Eventually, the AB distance will also be constant with time.

[[File:MT_DY.png|thumb|left|Momentum vs. time, dynamic calculation]]

[[File:IDT_DY2.png|thumb|center|Internuclear distance vs. time, dynamic calculation]] 

For the dynamic calculation, AB has an increased then constant momentum and BC has an oscillating momentum. This suggest that AB (ie. the speed of H1 moves away from H2-H3) increases at a constant speed once it reaches the 'flattened' potential surface and BC oscillates constantly. Thus, the AB distance becomes linearly proportional to time for large t values.

====Extra====
When the initial conditions and the final conditions are reversed, AB distance and BC distance are swapped, it is expected that the process will be exactly the reverse of the previous scenarios with the plots being the perfect mirror images of each other. 

[[File:Reversed_DT_yc16318.png|thumb|left|Internuclear Distance vs. time, Dynamic calculation, reversed]]

[[File:Reversed_MT_yc16318.png|thumb|center|Momenta vs. time, dynamic calculation, reversed]]

We can see from the plots that, while it resembles the reversed plot of the previous condition, it is not exactly the same. This could be due to the fact that the positions and the momenta of the final conditions were not perfectly estimated. This can also be supported by looking at its contour plot, where the trajectory doesn't pass the transition state due not enough energy in the system.

[[File:Reversed_SP_yc16318.png|thumb|center|Contour plot, dynamic calculation, reversed]]
 
activation energy calculation

===Q4===
Complete the table above 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? 
====Momenta vs. reactivity table====
{| class="wikitable" border="1"
|-
! p1/ g.mol-1.pm.fs-1 !! p2/ g.mol-1.pm.fs-1 !! Etot/ Kj.mol-1 !! Reactive? !! Description of the dynamics !! Illustration of the trajectory 
|-
| -2.56|| -5.1 || -414.280 || Yes || H3 approaches the non-vibrating H1-H2, forming a vibrating H2-H3 as H1 moves away. Products are formed. ||[[File:SP_C1-414_280.png|300px]]
|-
| -3.1 ||-4.1 || -420.077 || No || H3 collides with the noon-vibrating H1-H2. However, it does not move past the transition state as the potential energy exceeds the energy of the system. No products are formed, H3 moves away again. ||[[File:SP_C2-420_077.png|300px]]
|-
| -3.1 ||-5.1 || -413.977 || Yes || Similar to the first condition. H3 approaches the non-vibrating H1-H2, forming a vibrating H2-H3. The vibration is stronger than that in the first condition as the system possesses higher initial momentum. ||[[File:SP_C3-413_977.png|300px]]
|-
| -5.1 ||-10.1 || -357.277 || No || H3 approaches H1-H2. A H2-H3 is formed and vibrates strongly as H1 moves away. However, the system reverts back to the reactants and the H1-H2 is formed again. Eventually, H3 moves away, leaving a vibrating H1-H2 bond. The momentum of the system is large enough to allow the trajectory to cross the transition state twice, but no reaction happens in the end. The vibration is significantly large comparing to the previous conditions due to higher momentum. ||[[File:SP_C4-357_277.png|300px]]
|-
| -5.1 ||-10.6 || -349.477 || Yes ||H3 approaches H1-H2. At first, H2-H3 bond is formed and vibrates strongly as H1 moves away. However, the H1-H2 is not yet fully broken and vibrates once again. Eventually, a vibrating H2-H3 is formed as H1 moves away. The trajectory crosses the transition state three times before the products are formed. ||[[File:SP_C5-349_477.png|300px]]
|}

 

====Conclusion of the table====
From the table, we can conclude that Etot and the momenta of AB and BC are not the only factors for whether a reaction will happen. From the conditions in the table, we can see that, when p1 is kept constant, the condition with bigger p1/p2 difference is reactive, even though it has lower Etot. This could suggest that the difference between p1 and p2 might also play an important role in whether a system is reactive.

===Q5===
Given the results you have obtained, how will Transition State Theory predictions for reaction rate values compare with experimental values? 

In the transition state theory (TST), only the states of the reactants are considered for the rate of the reaction. It is assumed that, once the transition state is passed, the reaction will proceed to exit from the product channel irreversibly. Therefore, only one crossing of the transition state is allowed. Thus, the rate of the reaction can be approximated by the rate of molecules passing through the transition state.

However, as we can see from the previous simulations, it is possible for the reaction trajectories to pass the transition state twice or more if the particles possess large enough momenta. Therefore, we need a sufficiently large t to determine whether a certain pathway is reactive.

Additionally, because we are treating the nuclei as classical point masses, we expect the motions of the particle to be classical even though energy quantisation can be taken into account by the partition functions. Therefore, TST is unable to reflect the quantum effects in a system, such as the tunnelling effect.

London-Eyring-Polanyi-Sato (LEPS) potential 

== Exercise II: F - H - H system ==
===Q6===
By inspecting the potential energy surfaces, classify the F + H2 and H + HF reactions according to their energetics (endothermic or exothermic). How does this relate to the bond strength of the chemical species involved?
answer

===Q7===
Locate the approximate position of the transition state.
answer

===Q8===
Report the activation energy for both reactions.
answer

===Q9===
In light of the fact that energy is conserved, discuss the mechanism of release of the reaction energy. Explain how this could be confirmed experimentally.
answer

===Q10===
Discuss how the distribution of energy between different modes (translation and vibration) affect the efficiency of the reaction, and how this is influenced by the position of the transition state.
answer

File:Reversed SP yc16318.png

2020-05-08T16:12:13Z

Yc16318:

File:Reversed MT yc16318.png

2020-05-08T16:07:56Z

Yc16318:

2020-05-08T08:53:44Z

Yc16318: /* Momenta vs. reactivity table */