Abstract

The reaction profile (energy profile) is a widely used conceptual tool in chemical kinetics to represent the progress of a chemical reaction. Quantitatively, a reaction profile can be viewed as a minimum energy path (MEP) on the potential energy surface (PES), which connects the reactants and products through one or more transition states or intermediates. In this article, we used Mathematica program to demonstrate a generic method for finding reaction profile on a Müller-Brown PES by applying steepest descent algorithm. The properties of the MEP and stationary points were discussed in detail. The general characteristics of the transition state (TS), and imaginary mode were illustrated with a vibrational analysis of hydrogen exchange reaction, H₂+H → H+H₂.

1. Introduction

To accurately describe a reaction or its mechanism, an in-depth knowledge of PES of the system is very necessary. The PES gives us the energy of the system at a particular geometry or configuration, and other important derived quantities such as force (first derivative of energy with respect to the position, ∂E/∂r) or force constants (∂²E/∂r²) etc. In PES, the energy of the system (generally, the total energy), is represented as a function of a coordinate, (usually known as a reaction coordinate). For example one can find out the mechanism of a reaction (i.e. how and where the bond breaking and formation occur) by using molecular dynamics simulation in conjunction with the PES of the system. However, the construction of precise and detailed potential energy surfaces for a number of chemical reactions is not yet feasible. For example, consider a hydrogen exchange reaction. The PES of this reaction is well described in many standard text books in physical chemistry, and this PES is usually created by using hundreds of single-point ab initio calculations. The Figure 1, shows PES as a contour diagram of two independent variable R_ab and R_bc, and it runs between [0.0 - 3.0] Å. Here, a medium accurate step size of 0.10 Å was used, and we performed 900 (or 30²) single-point calculations to construct this PES (the curves in the figure can be drawn, for example, by using cublic-spline method).

And note that the H₂+H system is a simple system, and if we consider, say, H₂O+H reaction, (a bit more advanced system with twelve Cartesian degrees of freedom, DOF) the number of single point calculations to represent the PES in the Cartesian coordinate system with a step size of 0.10 Å would be 100¹², if we use a [0.0-10.0 Å] scale for each DOF. This example shows that one cannot obtain a complete PES of the system if the DOF is even twelve. However, there are other methods - which explore the chemically relevant portion of the PES ^{1, 2} exploring the evolution of the reactants to the products through the TS or intermediates along a least-action path. Applications of one of such schemes (known as intrinsic reaction coordinate (IRC)) ¹ using standard ab initio packages to some interesting reactions have been discussed erstwhile ^{3, 4, 5, 6, 7}. Note that IRC is just an another MEP, but defined in mass weighted coordinate system (See Appendix 1 in the supplementary information). For its applications in computational chemistry, refer Jensen’s text book ⁹.

In this article, we discuss 'minimum energy path' (MEP) algorithm ² in detail, since it represents a unique least-action path from the reactants to the products through the saddle points. Mathematically, MEP is defined as a steepest descent (SD) ⁸ path from the transition state (or a saddle point) towards a stationary point (equivalent to product, reactants or intermediates state since these points mathematically correspond to local or global minima on a surface). This SD can also be regarded as the energetically favorable (least-action) path between two minima.

Here, we demonstrate a generic MEP algorithm using a test (Müller-Brown) potential energy surface using Mathematica 10 program. Since we used only basic as well as common commands, one can also use equivalent computer algebra system such as Matlab, Maple, GNU-Octave etc. for doing this exercise by simply translating these commands. For brevity, we didn’t consider molecular dynamics calculations using this PES and we strongly suggest the reader to refer Sathyamurthy’s paper ¹⁰ for an introduction of MD using PES and other primitive concepts in the reaction dynamics.

The main steps in this exercise are following: (1) Construction of the PES, (2) Identifying stationary points (initial, intermediates, and final states) and transition states (TSs) etc., (3) Finding the imaginary mode at TS (4) Performing the MEP algorithm from TSs to obtain a complete reaction profile with an appropriate reaction coordinate. The stationary points and transition states were obtained by analytical differentiation. The ‘imaginary mode’ has been created by diagonalizing the force constant matrix (i.e. Hessian matrix) at the transition states. Forward and backward SD paths were calculated from the transition states to obtain the complete minimum energy path - which starts from the initial reactant location to the location of the final product via the saddle points and an intermediate state. Appendices are given finally, to discuss the properties of the Hessian matrix (Appendix 2 in the S.I.). Note that, a reader can download more detailed Mathematica 10 notebook or its PDF version from the links as shown in the supplementary information section.

2. Computational Procedures

Here we define MB-PES and its parameters (See Table 1). This PES is a well known test system used in Physical Chemistry for testing various numerical algorithms. The function, E(x,y), is stored in MBpes[x,y] variable and its plot is shown in the Figure 2. The closed (i.e. analytical) form of the MB-potential is also calculated. It should be noted that the MB-PES has two transition states and one intermediate state between the initial and final states (See below).

Now, we are about to explore the PES in detail. Arbitrarily, we can define that our initial state structure is approximately located at (-0.56x, 1.44y) and that of the product state structure is placed at (0.62x,0.02y). See the contour diagram in Figure 3 for the details. Similarly, one can easily locate two TSs and one intermediate state on the PES and these are located around the points at (-0.82x, 0.62y), (0.21x, 0.29y), and (-0.00x, 0.46y), respectively.

In order to calculate the MEP, we should find the exact figures of these points. Clearly, at these points we can write a simultaneous equation, viz.,

By applying the FindRoot function to these equations and from the approximate locations, we can calculate the exact numerical values of these points (the solutions are stored into some variables for later use, such as, xt, xbt, etc.). The exact position of TSs, initial, intermediate and final states are given in the Table 2.

All of these points can be now represented over the MB-PES contour diagram (See Figure 3).

To find the direction of descent from the TS to initiate the SD procedure, we need to find a starting vector and it is usually known as the transition vector or TV (also known as the imaginary mode, since it corresponds to a negative eigenvalue (See Appendix-2 for more details). This vector is obtained by diagonalizing the numerical Hessian matrix at the TS. For example, at TS1, the numerical Hessian matrix can be evaluated and its eigenvalues and eigenvectors are given in Table 3.

The above 2×2 matrix gives a single negative eigenvalue (or strictly, an imaginary frequency proportional to square root of -750.863) and that gives an imaginary mode, or eigenvector of [-0.761396 x + 0.648287 y]. Note that -1*[-0.761396 x +0.648287 y], is also an eigenvector of the Hessian matrix (and hence another TV). The above procedure of diagnolization can be applied to the second-Transition structure (ie. at (xbt, ybt)) and will give TVs at this point, viz., ±1*[-0.500306 x + 0.865849 y]. Table 4 contains the required the commands.

The resultant TV vectors at the TSs are shown below.

3. Running the MEP Algorithm from TS

The general equation for the steepest descent (SD) algorithm ³ is given as,

Here, R_i stands for a given position on SD path after a discrete SD step of i, ∇(R_i) is the gradient of the R_i and h is the step length (usually a very small number). Note that 'gradient of the R_i ' in our context is the force at R_i (i.e. ∂E/∂R). Starting from the transition vector and the TS geometry, one can trace the SD path towards adjacent minima. The TV is used only once (i.e. for the first step or when i =1). The rest of the SD steps use the gradient as shown in the above equation. Here we use the SD method to find the MEP towards the reactant. The result of the path is then stored in Pxb and Pyb arrays. For the sake of simplicity, we only demonstrate MEP from TS1 to the initial state and to the intermediate state (See Table 5 and Figure 4). And by extending these commands (which include the transition vector information at TS2), one can obtain the complete MEP of the reaction. See the Mathematica Notebook from the supplementary information for more details.

4. Results and Discussions

To get a meaningful MEP one should use an appropriate reaction coordinate - which represents the progress of the reaction. Since our PES is a function of Cartesian coordinates, an arbitrary function, f(x,y) can be used as a reaction coordinate and for simplicity, we chose a reaction coordinate that originates at the TS and is given by,

where r_i is the [x, y] coordinate of the MEP after the i'th SD step, and Δr_i is the Cartesian distance vector (Euclidian norm) from the two consecutive points on the path. This r_i is just the distance traveled along the MEP from the reference position to a particular point on the path.

By choosing TS1 as a reference point, we can divide the reaction profile into two parts; (1) MEP towards the reactant and (2) towards the product which passes through the intermediate and TS2. The following code (See, Table 6) calculates the 1-D reaction profile from the TS1 reference point.

This generic procedure can be used to determine, MEP, stationary points and TSs of analytical PES such as LEPS ¹² or other user generated PES with ab initio calculations. This algorithm can also be easily extended to higher dimensions by little modification of the scripts. The properties of the main locations are summarized in Table 7 (note that here r corresponds to the normal mode vectors). Finally, we would like to point that (a), for more accurate potential barrier one should consider zero point energy correction and, independently, the tunneling correction ¹³, and (b), estimation of MEP is widely used not only in gas phase molecular systems but also in higher dimensional hetero systems, eg. gas adsorption on metallic surfaces; See Refs. ¹⁴ and ¹⁵ for a review of some of the well accepted MEP methods that use more numerically efficient optimizing algorithms than SD, such as quasi-Newton ⁸ in this context. And in many of these methods, MEP calculation is initiated by using only initial and final state structures and not with the TS structure (as expected, locating TS is rather a difficult process in large dimensional systems). We strongly suggest the works of W.Quapp (Ref. ¹⁶ and references therein) for further reading on IRC and MEP and its algorithmics. A student can easily apply this analysis to LEPS potential ¹⁷ to find of the MEP of hydrogen exchange reaction or to ammonia inversion system ³. An insightful, albeit, mathematical description of 'least action path' is included in the classic book, Feynman lectures on Physics ¹⁸.

5. Discussions of MB-MEP

MEP gives a general information on the nature of reaction and mechanism. For example, concerning the calculated energy profile one can easily understand that the overall reaction is endothermic. This energy landscape gives another valuable information that the reaction proceed via forming an intermediate state (in practical situation its life span would be a transient or very small). Hammond’s postulates can be applied here ¹¹ and it indicates the characteristics of the geometries at TSs (for example, TS1 geometry will be similar to that of the initial structure). Since TS1 looks much more like reactants than the final product, it can be considered as an early barrier system, so a translation energy (heat) can be used, instead of vibrational excitation, to promote the reaction (Polanyi’s rule, ¹²).

Moreover, one can easily identify activation energy (or barrier, E₁), the enthalpy of the reaction (ie. enthalpy of formation of the product, E₃), and E₂ - the activation barrier for the formation of the intermediate. They are all are represented in the Figure 6.

Using this energy profile, we can easily find the reaction rate constant, k, by using Arrhenius Equation, k = A e^-E₁^RT (See ¹² for more discussions). Here, A is the pre-exponential factor, R and T are the universal gas constant and temperature in Kelvin, respectively. Note that, A, cannot be determined from the minimum energy path alone. The most widely used approach to finding A is transition state theory, which estimates A on the basis of the shapes of the PES in the vicinities of the reactant and the transition state for each step. If the PES near the transition state is broad and flat, transition state theory will give a large A, while if it is narrow and has steep walls, the theory will yield a small A. Additionally, one can analysis the evolution of geometries (as well as with molecular orbitals) of the reactants along this profile and it will give an idea of the nature of reaction mechanism, in detail.

The transition vector - the direction of the SD path at TS – represents the reaction coordinate which give the progress of the reaction. To illustrate this see the below figure, which depicts the vibrational modes of the TS of the exchange reaction of the H₂+H.

Note that, the TS is a 3 atom, linear molecular system. So, there are (3_˟3-5) modes are possible (i.e. 4 modes). We have used Gaussian STO-3G calculation to find the optimized TS and performed vibrational calculations for this geometry. And it clearly gives 4 frequencies with one imaginary frequency (imaginary frequency is usually represented by a negative frequency, ω = -3026cm^-1, or more correctly 3026i cm^-1; since ω = where k, and m are the force constant and effective mass of the system, respectively). Here, the force constant (actually it is a 1_˟1 force constant matrix) for this particular vibration can be represented by:

and this curvature of the PES obviously represents negative parabola (see the first curve in the Figure 7) which passes through the TS and it makes the MEP. The positions, P1 and P2 of the parabola, will be given by the transition vector (TV for P1 and –TV for the P2, for example). All the three other parabolas of the other three vibrations are the positive parabolas in the PES. Note that the TV clearly indicates the direction of the reaction (bond dissociation and formation). Furthermore, steepness of the parabola represents that magnitude of the vibration, in other words, parabola with a greater steepness will results stronger vibrations.

6. Conclusion

In this Mathematica demonstration we illustrated the MEP and its reaction profile by using MB-PES analytical PES. Stationary points and TSs were located and from the transition vectors at TS1 and TS2, MEP was obtained by applying steepest descent algorithm. The MEP summarizes information about the reaction mechanism. The presence of an intermediate is very clear on the reaction profile diagram. Because the potential energy described by the MEP is usually the most important contributing factor to the internal energy or enthalpy of the reacting system, it is a reasonable approximation to relate potential energy differences along the MEP to changes of energy or enthalpy of the chemical system. Reaction coordinate was calculated and it was used to plot 1D reaction profile. Finally the transition state structure of hydrogen exchange reaction was analyzed and the meaning of the imaginary mode was discussed in terms of the reaction coordinate.

References

Supplementary Information

Mathematica (v10) Notebook of this article is available at WJCE web site.
SupportingInformation-MEP_on_MBsurface-Mathematica10.pdf
SupportingInformation-MEP_on_MBsurface-Mathematica10.nb