World Journal of Chemical Education
Volume 10, 2022 - Issue 4
Website: https://www.sciepub.com/journal/wjce

ISSN(Print): 2375-1665
ISSN(Online): 2375-1657

Article Versions

Export Article

Cite this article

- Normal Style
- MLA Style
- APA Style
- Chicago Style

Research Article

Open Access Peer-reviewed

Krishnamohan G. P.^{ }, Omar H., Sreeja T. D., Roy K. B.

Received July 26, 2022; Revised September 04, 2022; Accepted September 12, 2022

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_{2}+H → H+H_{2}.

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 (∂^{2}E/∂**r**^{2}) 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

**Figure****1****.**PES of the hydrogen exchange reaction obtained from a set of ab initio calculations; the asterisk indicates the TS and unit of bond length is in Å

And note that the H_{2}+H system is a simple system, and if we consider, say, H_{2}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^{12}, 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)) ^{ 1} 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 ^{ 9}.

In this article, we discuss *'minimum energy path'* (MEP) algorithm ^{ 2} 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) ^{ 8} 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 ^{ 10} 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.

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).

**Figure****2****.**The above figure shows the TSs, intermediate (I), and stationary points (initial, S_{I}and final, S_{F}) on the generated MB-PES

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

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.

The general equation for the *steepest descent (SD) *algorithm ^{ 3} 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.

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

where **r**_{i} is the [x, y] coordinate of the MEP after the * i*'th SD step, and Δ

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 ^{ 12} 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 (

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 ^{ 11} 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, ^{ 12}).

Moreover, one can easily identify activation energy (or barrier, E_{1}), the enthalpy of the reaction (ie. enthalpy of formation of the product, E_{3}), and E_{2} - 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}_{1}^{RT} (See ^{ 12} for more discussions). Here, *A** *is the pre-exponential factor,

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_{2}+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 3026*i* 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.

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.

[1] | Fukui, K, Nobel lecture, https://www.nobelprize.org/uploads/ 2018/ 06/fukui-lecture.pdf. | ||

In article | |||

[2] | Henkelman, G. and Jonsson, H, “Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points”, J. Chem. Phys., 113, 9978-9985, November 2000. | ||

In article | View Article | ||

[3] | Halpern, A. M., Ramachandran, B. R. and Glendening, E, “The Inversion Potential of Ammonia: An Intrinsic Reaction Coordinate Calculation for Student Investigation”, J. Chem. Educ., 84, 1067-1072, June 2007. | ||

In article | |||

[4] | Goldstein, E. and Leach, A. G, “Energy Contour Plots: Slices through the Potential Energy Surface That Simplify Quantum Mechanical Studies of Reacting Systems”, J. Chem. Educ., 83, 451-459, March 2006. | ||

In article | View Article | ||

[5] | Halpern, A. M, “Computational Studies of Chemical Reactions: The HNC–HCN and CH_{3}NC–CH_{3}CN Isomerizations”, J. Chem. Educ., 83, 69-76, January 2006. | ||

In article | |||

[6] | Lehman, J. J. and Goldstein, E, “The Potential Energy Surface of ClF_{3}”, J. Chem. Educ., 73, 1096-1098, November 1996. | ||

In article | |||

[7] | Marzzacco, C. J.and Baum, J. C, “Computational Chemistry Studies on the Carbene Hydroxymethylene”, J. Chem. Educ., 88, 1667-1671, August 2011. | ||

In article | View Article | ||

[8] | Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P, “ Numerical Recipes in C: Art of Scientific Computing”, Cambridge University Press, 1999. | ||

In article | |||

[9] | Jensen F, “Introduction to Computational Chemistry”, 3rd Edition, Wiley, 2017. | ||

In article | |||

[10] | Sathyamurthy, N. and Joseph, T, “Potential energy surface and molecular reaction dynamics”, J. Chem. Educ., 61, 968-971, November 1984. | ||

In article | |||

[11] | Carey, F. A., and Sundberg, R .J, “Advanced Organic Chemistry, Part A: Structure and Mechanisms”, Springer; 4th edition, 2006. | ||

In article | |||

[12] | Laidler, K. J, “Chemical Kinetics”, Pearson; 3rd edition, 1997. | ||

In article | |||

[13] | Truhlar, D (Ed.), “Potential Energy Surfaces and Dynamics Calculations: for Chemical Reactions and Molecular Energy Transfer”, Springer, 2013. | ||

In article | |||

[14] | Sheppard, D., Terrell, R. and Henkelman, G, “Optimization methods for finding minimum energy path”, J. Chem. Phys., 128, 134106, April 2008. | ||

In article | View Article PubMed | ||

[15] | Henkelman, G., Jóhannesson, G. and Jónsson, H, “Theoretical Methods in Condensed Phase Chemistry. Progress in Theoretical Chemistry and Physics”, vol. 5, Springer, Dordrecht, 2002. | ||

In article | |||

[16] | Quapp, W, “Chemical Reaction Paths and Calculus of Variations”, Theo. Chem. Acc., 121, 227–237, July 2008. | ||

In article | |||

[17] | Moss, S. J. and Coady, C. J, “Potential energy surfaces and transition state theory”, J. Chem. Educ., 60, 455, June 1983. | ||

In article | View Article | ||

[18] | Feynman, R, “Feynman lectures on Physics”, vol. 2, sec.19, Pearson, 2012. | ||

In article | |||

Mathematica (v10) Notebook of this article is available at WJCE web site.

SupportingInformation-MEP_on_MBsurface-Mathematica10.pdf

SupportingInformation-MEP_on_MBsurface-Mathematica10.nb

Published with license by Science and Education Publishing, Copyright © 2022 Krishnamohan G. P., Omar H., Sreeja T. D. and Roy K. B.

This work is licensed under a Creative Commons Attribution 4.0 International License. To view a copy of this license, visit https://creativecommons.org/licenses/by/4.0/

Krishnamohan G. P., Omar H., Sreeja T. D., Roy K. B.. Exploring Potential Energy Surface with *Mathematica*: An Algorithmic Demonstration of Minimum Energy Path, Stationary Points and Transition State. *World Journal of Chemical Education*. Vol. 10, No. 4, 2022, pp 124-130. https://pubs.sciepub.com/wjce/10/4/1

P., Krishnamohan G., et al. "Exploring Potential Energy Surface with *Mathematica*: An Algorithmic Demonstration of Minimum Energy Path, Stationary Points and Transition State." *World Journal of Chemical Education* 10.4 (2022): 124-130.

P., K. G. , H., O. , D., S. T. , & B., R. K. (2022). Exploring Potential Energy Surface with *Mathematica*: An Algorithmic Demonstration of Minimum Energy Path, Stationary Points and Transition State. *World Journal of Chemical Education*, *10*(4), 124-130.

P., Krishnamohan G., Omar H., Sreeja T. D., and Roy K. B.. "Exploring Potential Energy Surface with *Mathematica*: An Algorithmic Demonstration of Minimum Energy Path, Stationary Points and Transition State." *World Journal of Chemical Education* 10, no. 4 (2022): 124-130.

Share

[1] | Fukui, K, Nobel lecture, https://www.nobelprize.org/uploads/ 2018/ 06/fukui-lecture.pdf. | ||

In article | |||

[2] | Henkelman, G. and Jonsson, H, “Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points”, J. Chem. Phys., 113, 9978-9985, November 2000. | ||

In article | View Article | ||

[3] | Halpern, A. M., Ramachandran, B. R. and Glendening, E, “The Inversion Potential of Ammonia: An Intrinsic Reaction Coordinate Calculation for Student Investigation”, J. Chem. Educ., 84, 1067-1072, June 2007. | ||

In article | |||

[4] | Goldstein, E. and Leach, A. G, “Energy Contour Plots: Slices through the Potential Energy Surface That Simplify Quantum Mechanical Studies of Reacting Systems”, J. Chem. Educ., 83, 451-459, March 2006. | ||

In article | View Article | ||

[5] | Halpern, A. M, “Computational Studies of Chemical Reactions: The HNC–HCN and CH_{3}NC–CH_{3}CN Isomerizations”, J. Chem. Educ., 83, 69-76, January 2006. | ||

In article | |||

[6] | Lehman, J. J. and Goldstein, E, “The Potential Energy Surface of ClF_{3}”, J. Chem. Educ., 73, 1096-1098, November 1996. | ||

In article | |||

[7] | Marzzacco, C. J.and Baum, J. C, “Computational Chemistry Studies on the Carbene Hydroxymethylene”, J. Chem. Educ., 88, 1667-1671, August 2011. | ||

In article | View Article | ||

[8] | Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P, “ Numerical Recipes in C: Art of Scientific Computing”, Cambridge University Press, 1999. | ||

In article | |||

[9] | Jensen F, “Introduction to Computational Chemistry”, 3rd Edition, Wiley, 2017. | ||

In article | |||

[10] | Sathyamurthy, N. and Joseph, T, “Potential energy surface and molecular reaction dynamics”, J. Chem. Educ., 61, 968-971, November 1984. | ||

In article | |||

[11] | Carey, F. A., and Sundberg, R .J, “Advanced Organic Chemistry, Part A: Structure and Mechanisms”, Springer; 4th edition, 2006. | ||

In article | |||

[12] | Laidler, K. J, “Chemical Kinetics”, Pearson; 3rd edition, 1997. | ||

In article | |||

[13] | Truhlar, D (Ed.), “Potential Energy Surfaces and Dynamics Calculations: for Chemical Reactions and Molecular Energy Transfer”, Springer, 2013. | ||

In article | |||

[14] | Sheppard, D., Terrell, R. and Henkelman, G, “Optimization methods for finding minimum energy path”, J. Chem. Phys., 128, 134106, April 2008. | ||

In article | View Article PubMed | ||

[15] | Henkelman, G., Jóhannesson, G. and Jónsson, H, “Theoretical Methods in Condensed Phase Chemistry. Progress in Theoretical Chemistry and Physics”, vol. 5, Springer, Dordrecht, 2002. | ||

In article | |||

[16] | Quapp, W, “Chemical Reaction Paths and Calculus of Variations”, Theo. Chem. Acc., 121, 227–237, July 2008. | ||

In article | |||

[17] | Moss, S. J. and Coady, C. J, “Potential energy surfaces and transition state theory”, J. Chem. Educ., 60, 455, June 1983. | ||

In article | View Article | ||

[18] | Feynman, R, “Feynman lectures on Physics”, vol. 2, sec.19, Pearson, 2012. | ||

In article | |||