MRD:01507626

H+H₂ system

The initial conditions of the simulation are shown in the table below.

	Distance/pm	Momentum/g.mol-1.pm.fs-1
AB	230	-5.2
BC	74	0

Figure 1. Potential energy surface and trajectory of H+H2 reaction

The results of the calculation are shown in Figure 1 in the form of the potential energy surface. The black line at the bottom of the surface shows the trajectory taken. The line is not straight due to the osciallations of the diatomic molecules. The part of the trajectory in which the BC distance stays constant at around 74 pm and AB decreases to about 74 pm, represents the reactants. The corner in the trajectory is the transition state of the reaction. The part where AB stays constant at 74pm while BC increases shows atom C moving away after the reaction while BC remains.

Transition state

The transition state is a saddle point on the potential energy surface (PES). Saddle points, maxima and minima are all stationary points. For a point to be a stationary point of the potential energy surface all partial derivatives have to equal 0 at the point.

${dV \over dr_i} = 0$

The nature of the stationary point can then be determined by expanding the n-dimensional PES using the Taylor series around the point where N is the number of atoms. In most cases the quadratic approximation described below is suitable.

By introducing the Hessian matrix M one can rewrite the expansion.

The Hessian matrix is real and symmetric and therefore has n real eigenvalues λ_k and n orthogonal eigenvectors e_k.

ẟ_rs is Kronecker delta, showing the orthogonality between the eigenvectors. Since one can construct an n-dimesional vecor out of the n eigenvectors as base vectors. Therefore Δr can be expanded in terms of the eigenvectors.

Substituing the last equation into the taylor expansion, one obtains the following. Another way of looking at this substitution is decomposing the vector Δr into the eigenvectors of M.

One can think of the result as determining the nature of the stationary point along one eigenvector of M, along one dimension in the n dimensional vector speace. A positive eigenvalue represents a minimum in a given dimension while a negative eigenvalue represents a maximum for a given dimension. This leads us to the conclusion that if all eigenvalues are positive for all coefficients a_k, the point is a minimum. If all eigenvalues are negative for all coefficients a_k the point is a maximum. If some eigenvalues are positive while others negative the point is a saddle point. In the case of a 2 variable function a saddle point is defined as one eigenvalue having a different sign to the other, meaning along one eigenvector of M it is a minimum while a maximum along another orthogonal to the first one. In the case of an n variable function the point can be a maximum along an arbitrary number of eigenvectors, k, where k smaller than n-1, while being a minimum along n-k other orthogonal eigenvectors. However, a transition state only corresponds to those points where one eigenvalue is negative while all others positive. This ensures that it is still a point along the minimum energy path.

It should be noted that this method may fail in case some eigenvalues are 0 while others have the same sign. In this case numerical approximate methods need to be used.

In practice, the full PES is rarely computed due to the large computational cost and algorithms like the Berny algorithm need to be used. It starts off with a well guessed starting geometry, computes the Hessian matrix and takes a step on the PES towards the transition state and lowering the energy gradient, the condition of a stationary point. This method is only feasible when a good starting geometry was provided. This can be obtained by scanning along one direction on the PES, for example the by freezing all bond legths except one.

Locating the transition state

When looking at figure 2 one can see that the PES is symmetric along the plane r_AB=r_BC. Therefore The transition state must have the same symmetry. One way to search for the transition state is to start with the condition r_AB=r_BC and vary the distance. At the transition state, the particles should not move as long as no initial momentum is given. In this specific system transition state search is greatly simplified due to the constraint given by symmetry. If the PES was known the method of Lagrange undetermined multipliers could be used to determine an analytic answer. However, as the PES is not fully known the easiest method would be to find the minimum along the constraint. There are several methods to doing this but the simplest is using gradient descent along the constraint. The formula of the steepest descent is the following where k is a pre defined step length.

$r_{n+1}-r_n=-k\frac{\nabla V(r_n)}{\lVert {\nabla V(r_n)}\rVert}$

Using the definition of force as the derivative of the energy this can be rewritten.

$r_{n+1}-r_n=-k\frac{F(r_n)}{\lVert {F(r_n)}\rVert}$

Figure 3. Internuclear distance vs time of initial trajectory

From figure 3 we can see that 90pm is a valid first guess. The step length was chosen as 0.5 pm. Convergence is reached when the sign of force vector changes. Then the step size is reduced to 0.05 pm.

n	rn	Fn
0	90	$\begin{pmatrix} 0.132 \\ 0.132 \end{pmatrix}$
1	90.3536	$\begin{pmatrix} 0.071 \\ 0.071 \end{pmatrix}$
2	90.7072	$\begin{pmatrix} 0.011\\ 0.011 \end{pmatrix}$
3	91.0608	$\begin{pmatrix} -0.047\\ -0.047 \end{pmatrix}$

The transition state lies between 90.7072 pm and 91.0608 pm. The starting point is chosen by linearly interpolating the distance as 90.7743. The forces at this point are computed as 0 in both directions. Therefore the transition state has been reached at 90.7743. By interpolating the steps before and after the sign change, the converges much faster. For higher dimensional optimisation problems gradient descent is not very efficient close to the stationary point but efficient at optimisation far from the stationary point. Methods like Newton-Raphson or the previously discussed Berny-Algorithm that follows eigenvectors are much more efficient in those cases. The separation of 90.7743 pm can be confirmed as a transition state with the help of figure 4. No oscillations in the atom separations show that the atoms are still at their positions.

Figure 4. Internuclear distance vs time plot when initial positions are set corresponding to the transition state.

Comparison of minimum energy trajectory (mep) with Dynamic trajectory

The mep corresponds to the reaction path with infinitely slow motion. The momentum of the atoms is set to 0 at every time step so that the trajectory follows the floor of the potential energy surface as shown in figure 5. When changing the type of the calculation back to dynamic as shown in figure 6, one can clearly see the vibrational energy of the products. In both cases the initial conditions were set as r_AB=90.7 pm and r_BC=90.7743 pm.

Trajectories near the transition state

	rAB=rBC=rTS=90.7743	rAB=91.7743 and rBC=rTS=90.7743	rAB=322.05 and rBC=74.34, pAB=-5.062 and pBC=-1.887
Internuclear distance vs Time		Final rAB=322.05 Final rBC=74.34	Final rAB=92.1 Final rBC=90.2
Momentum vs time		Final pAB=5.062 Final pBC=1.887	Final pAB=0.066 Final pBC=-0.072

One can see that at the transition state there is no displacement, the atoms are still at their positions. A slight displacement of 1 pm from the transition state causes the distance r_BC to decrease and oscillate while atom A steadily moves away from the molecule shown by the increasing distance r_AB. When using the final state of the system for the initial state of another simulation with the signs of the momenta reversed, One can see that it follows along the same trajectory as before in reverse order. The reason why the transition state is not perfectly reached may be due to the uncertainty with which the final positions and momenta were determined.

Reactive and unreactive trajectories

In the following table a number of trajectories are analysed. In all cases the original starting distances were set to r_AB=200 pm and r_BC=74 pm.

pBC/ g.mol-1.pm.fs-1	pAB/ g.mol-1.pm.fs-1	Etot	Reactive?	Description of the dynamics	Illustration of the trajectory
-2.56	-5.1	-414.280	Yes
-3.1	-4.1	-420.077	No
-3.1	-5.1	-413.977	Yes
-5.1	-10.1	-357.277	No
-5.1	-10.6	-349.477	Yes

Transition state theory

Transition state theory is based on three key assumptions.

1. All systems with kinetic energy in the direction of the reaction coordinate greater than activation energy will be reactive.
2. Kinetic energy along the reaction coordinate is populated according to the Boltzmann distribution.

When looking at the fourth trajectory in the previous section one can clearlz see that the system has crossed the transition state, entered the prodcut channel, crossed the transition state again and returned to the reactant channel. It is therefore an unreactive trajectory. However, according to Transition state theory it would be counted as a reactive one. Transition state theory therefore overestimates the rate of reaction. Transition state theory also doesn't take into account quantum tunneling since it is a classical theory. Quantum tunneling is an effect that increases the rate of reaction compared to transition state theory because a system can react by never reaching the transition state and tunneling through the activation barrier. For most systems this effect is however negligible since quantum is only substantially observed for very light particles. For some reactions this effect can be considerable but in most cases it will be negligible. Therefore, generally speaking transition state theory overestimates the rate of reaction.

F+H₂ system

Inspection of the PES

The following two figures show the trajectories of the F+H₂ and the HF+F reaction respectively.

One can clearly see that for the F+H₂ reaction the reactant channel is much higher than the product channel with an early transition state. One can clearly see the large vibrational excitation in the product channel. In the HF+H reaction the reactant channel is much lower than the product channel. The H atom has to collide with the molecule with a large momentum to provide the system with enough energy to overcome the transition state. The vibrations of the resulting H₂ molecule are much smaller as seen in the Figure on the right.

Locating of the transition state

Hammond's postulate states that the transition state resembles the reactants or products of a reaction depending on whichever one it is closest to in energy. When selecting the exothermic F+H₂ reaction from Hammond's postulate one can deduce that the transition state should resemble the reactants. When looking at the contours of the potential energy surface one can see that this is in fact the case. lower energy state of the HF+H state compared to the F+H₂ state indicates that the HF bond is much stronger than the H₂ bond.

When it comes to finding the transition state precisely, various methods exist. Previously, the gradient descent method could be used due to the high symmetry of the system. In this case a two dimensional search method needs to be used. When faced with such a problem for a larger system, the initial guess for the transition state is obtained by the use of synchronous transit methods that estimate the position of the transition state either linearly or quadratically along a straight line or curve between reactants and products. Further search is then carried out using other methods like adapted Newton-Raphson or eigenvector following. In the following the latter algorithm will be employed to find the transition state. The basic principle of the algorithm comes from an intuition about the eigenvectors of the Hessian matrix introduced earlier. The eigenvectors point into a direction in which the curvature is independent of any other directions. The sign of the eigenvalue corresponds to the sign of the curvature in the direction of the eigenvector. From the definition of a transition state given earlier one could find a transition state by minimizing one eigenvalue while maximizing all others. The first step of the RFO algorithm is to determine the shift parameter λ. In the simple gradient descent parameter that was previously used to obtain the transition state position this parameter was simply set to a sensible value. After rewriting the information included in the quadratic Taylor approximation in the form of a rational function with the help of a scaling matrix S. S is taken as the identity in the following since the scaling can be absorbed into the coordinates.

$\Delta V=\frac{g^T\Delta r+0.5\Delta r^T M \Delta r}{1+\Delta r^T S \Delta r}$

$\begin{pmatrix} M & \nabla \Delta r \\ (\nabla (\Delta r)^T & 0 \end{pmatrix}\begin{pmatrix} \Delta r\\ 1 \end{pmatrix}=\lambda \begin{pmatrix} S & 0 \\ 0 & 1 \end{pmatrix}\begin{pmatrix} (\Delta r)\\ 1 \end{pmatrix}$ where λ=2Δr

One can represent this equation as two separate equations.

$(M-\lambda S)\Delta r+g=0$

$g^T\Delta r=\lambda$ where g is the gradient of the PES at r_i

The equation can be rewritten in component form.

$\Delta r_k=\frac{g_k}{\lambda - m_k}$ where m_k are the eigenvalues of M and the diagonal elements of the diagonalised Hessian M

We can substitute this equation back into the separated RFO equation and obtain an expression for the change in V.

$\Delta V=\frac{1}{\sum_{k}r_k} \sum_{k}g_k \frac{\lambda - 0.5m_k}{{(\lambda - m_k)}^2}$

Therefore a step along or against the gradient in the direction of a specific eigenvector is determined by the sign of λ - m_k. One possible way of using this algorithm is to determine the lowest eigenvalue of the original RFO matrix, compute the gradient and use the formula for Δr_i. Another method would be to use the second lowest eigenvalue of M as λ. Therefore the algorithm takes steps upwards in the direction of one eigenvector while simultaneously taking steps downwards in the direction of all other eigenvectors. In this form eigenvector following can be used to find a first order saddle point.

As initial initial guess for the search of the transition state on the F-H-H PES r_FH=190 and r_HH=75.

Step	Force, g	λ	m1,m2	rFH,rHH
0	$\begin{pmatrix} 0.006\\ -0.224 \end{pmatrix}$	0.323	-0.001,0.323

Reaction dynamics

F+H₂ reaction

Initial conditions	Momentum vs time	Surface plot	Description of Dynamics
rHH = 74 pm rFH = 200 pm, pHH = 0 g.mol-1.pm.fs-1 pFH = -1.0 g.mol-1.pm.fs-1
rHH = 74 pm rFH = 200 pm, pHH = -6.1 g.mol-1.pm.fs-1 pFH = -1.0 g.mol-1.pm.fs-1
rHH = 74 pm rFH = 200 pm, pHH = -5.1 g.mol-1.pm.fs-1 pFH = -1.0 g.mol-1.pm.fs-1
rHH = 74 pm rFH = 200 pm, pHH = -4.1 g.mol-1.pm.fs-1 pFH = -1.0 g.mol-1.pm.fs-1
rHH = 74 pm rFH = 200 pm, pHH = -3.1 g.mol-1.pm.fs-1 pFH = -1.0 g.mol-1.pm.fs-1
rHH = 74 pm rFH = 200 pm, pHH = -2.1 g.mol-1.pm.fs-1 pFH = -1.0 g.mol-1.pm.fs-1
rHH = 74 pm rFH = 200 pm, pHH = -1.1 g.mol-1.pm.fs-1 pFH = -1.0 g.mol-1.pm.fs-1
rHH = 74 pm rFH = 200 pm, pHH = 1.1 g.mol-1.pm.fs-1 pFH = -1.0 g.mol-1.pm.fs-1
rHH = 74 pm rFH = 200 pm, pHH = 2.1 g.mol-1.pm.fs-1 pFH = -1.0 g.mol-1.pm.fs-1
rHH = 74 pm rFH = 200 pm, pHH = 3.1 g.mol-1.pm.fs-1 pFH = -1.0 g.mol-1.pm.fs-1
rHH = 74 pm rFH = 200 pm, pHH = 4.1 g.mol-1.pm.fs-1 pFH = -1.0 g.mol-1.pm.fs-1
rHH = 74 pm rFH = 200 pm, pHH = 5.1 g.mol-1.pm.fs-1 pFH = -1.0 g.mol-1.pm.fs-1
rHH = 74 pm rFH = 200 pm, pHH = 6.1 g.mol-1.pm.fs-1 pFH = -1.0 g.mol-1.pm.fs-1
rHH = 74 pm rFH = 200 pm, pHH = 0.2 g.mol-1.pm.fs-1 pFH = -1.6 g.mol-1.pm.fs-1

HF+H reaction

Initial conditions	Momentum vs time	Surface plot	Description of Dynamics
rHH = 200 pm rFH = 91 pm, pHH = -17 g.mol-1.pm.fs-1 pFH = 0 g.mol-1.pm.fs-1
rHH = 200 pm rFH = 91 pm, pHH = -17 g.mol-1.pm.fs-1 pFH = 0.1 g.mol-1.pm.fs-1