Dear Label_Name,

Greetings! I hope you are well.

I am writing on behalf of an international publishing house, Science and Education Publishing, USA (http://www.sciepub.com).

Science and Education Publishing (SciEP) is a multidisciplinary open access publisher of Journals covering the fields of Science, Technology and Medicine. All SciEP publications are subject to high-quality peer review, editorial and production processes and are freely available online.

Welcome to the April, 2013 edition of SciEP newsletter, providing monthly summaries of articles on science & education related topics of interest to scientist and professionals. The mission of our journals is to provide an open access forum for publishing new findings on issues pertinent to scientific and academic field. We also invite you to submit original research papers, review articles, short reports, or case studies are considered for publication in our Multidisciplinary Journals.

American Journal of Public Health Research
DOI: 10.12691/ajphr-1-1-3

The Implementation of a Tele-Homecare System with Service Oriented Architecture and HL7 Message Transmission Standard
Ching-Sung Wang

American Journal of Medical and Biological Research
DOI: 10.12691/ajmbr-1-1-2

The Influence of Ethanol on Pyruvate Kinases Activity in Vivo, in Vitro, in Silico
Lelevich Sergey Vladimirovich, Khrustalev Vladislav Victorovich, Barkovsky Eugene Victorovich, Shedogubova Tatyana Aleksandrovna

American Journal of Medical Sciences and Medicine
DOI: 10.12691/ajmsm-1-1-1

Clinical and Histological Features in Wegener Granulomatosis
Li Cavoli Gioacchino, Passantino Rita, Bono Luisa, Tortorici Calogera, Ferrantelli Angelo, Giammarresi Carlo, Rotolo Ugo

American Journal of Infectious Diseases and Microbiology
DOI: 10.12691/ajidm-1-1-2

The Curse of Lymphatic Filariasis: Would the Continual Use of Diethylcarbamazine Eliminate this Scourge in Papua New Guinea?
Anne Lanham, Lillian Mwanri

We have released issues of some journals on February 2013.
American Journal of Public Health Research
American Journal of Medical Sciences and Medicine
American Journal of Medical and Biological Research
American Journal of Pharmacological Sciences
American Journal of Clinical Medicine Research
American Journal of Nursing Research
American Journal of Cancer Prevention
American Journal of Medicine Studies
American Journal of Infectious Diseases and Microbiology
American Journal of Cardiovascular Disease Research
American Journal of Biomedical Research
American Journal of Mechanical Engineering
American Journal of Microbiological Research
Journal of Preventive Medicine
American Journal of Food and Nutrition
American Journal of Hypertension Research
American Journal of Sports Science and Medicine
Journal of Cancer Research and Treatment
American Journal of Zoological Research
Journal of Biomedical Engineering and Technology
American Journal of Microbiological Research
Biomedical Science and Engineering

We are pleased to announce that submissions are now being accepted for peer review for all the journals across all our subject areas in SciEP.

Please submit your manuscript via our online submission system:
http://www.sciepub.com/author/signin

To unsubscribe from our mailing list, please click on the following link: unsubscribe

I am looking forward to hearing from you.

Kind regards,

David Rain
Email:join@sciepub.com
Science and Education Publishing (http://www.sciepub.com)

Exploring Potential Energy Surface with Mathematica: An Algorithmic Demonstration of Minimum Energy Path, Stationary Points and Transition State

Krishnamohan G. P.*, Omar H., Sreeja T. D., Roy K. B.

Department of Science and Humanities, Mar Baselios College of Engineering and Technology (autonomous), Nalanchira, Trivandrum, Kerala 695015, India

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

Keywords: minimum energy path, potential energy surface, transition state, stationary state, imaginary mode

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

Citation: 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, 10 (4). 124-130. January 2022. doi:10.12691/wjce-10-4-1

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 Rab and Rbc, 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).

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₂+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)) [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.

2. Computational Procedures

1) Construction of the PES

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

Table 1. Construction of the MB-PES

Figure 2. The above figure shows the TSs, intermediate (I), and stationary points (initial, SI and final, SF) on the generated MB-PES

2) Identifying Stationary Points

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.

Table 2. Finding the exact positions of the TS and stationary points

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

3) Finding the Imaginary Modes at TSs

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.

Table 3. Imaginary mode of TS1

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.

Table 4. Imaginary mode of TS2

The resultant TV vectors at the TSs are shown below.

Figure 3. MB-PES contour diagram (Blue, Cyan, Red and Green dots represents, the initial, intermediate, TSs and final points of the reaction, respectively) with transition vectors at TS1 and TS2 (red arrows)

3. Running the MEP Algorithm from TS

The general equation for the steepest descent (SD) algorithm [3] is given as,

Here, Ri stands for a given position on SD path after a discrete SD step of i, ∇(Ri) is the gradient of the Ri and h is the step length (usually a very small number). Note that 'gradient of the Ri ' in our context is the force at Ri (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.

Table 5. MEP algorithm demonstration from the TS1 to reactant and to the intermediate

Figure 4. MEP from the TS1 to the reactant (first figure) and to the intermediate (second figure). The overall figure is shown in the third figure. The last figure is the complete MEP from the initial state to the final state

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 ri is the [x, y] coordinate of the MEP after the i'th SD step, and Δri is the Cartesian distance vector (Euclidian norm) from the two consecutive points on the path. This ri 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.

Table 6. Complete MEP from the TS1 state towards reactant and the product positions

Figure 5. The whole 1D reaction profile calculated from TS1 as a function of r(x,y)

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 (a), for more accurate potential barrier one should consider zero point energy correction and, independently, the tunneling correction [13], 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. [14] and [15] for a review of some of the well accepted MEP methods that use more numerically efficient optimizing algorithms than SD, such as quasi-Newton [8] 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. [16] and references therein) for further reading on IRC and MEP and its algorithmics. A student can easily apply this analysis to LEPS potential [17] to find of the MEP of hydrogen exchange reaction or to ammonia inversion system [3]. An insightful, albeit, mathematical description of 'least action path' is included in the classic book, Feynman lectures on Physics [18].

Table 7. Properties of the main locations at the PES

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 [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₁), 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.

Figure 6. Reaction profile with some important parameters on this

Using this energy profile, we can easily find the reaction rate constant, k, by using Arrhenius Equation, k = A e⁻E₁RT (See [12] 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⁻¹, or more correctly 3026i cm⁻¹; 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.

Figure 7. Four different vibrational modes of the TS of the hydrogen exchange reaction

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.

[1]   | Fukui, K, Nobel lecture, https://www.nobelprize.org/uploads/ 2018/ 06/fukui-lecture.pdf.
[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.
[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.
[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.
[5]   | Halpern, A. M, “Computational Studies of Chemical Reactions: The HNC–HCN and CH₃NC–CH₃CN Isomerizations”, J. Chem. Educ., 83, 69-76, January 2006.
[6]   | Lehman, J. J. and Goldstein, E, “The Potential Energy Surface of ClF₃”, J. Chem. Educ., 73, 1096-1098, November 1996.
[7]   | Marzzacco, C. J.and Baum, J. C, “Computational Chemistry Studies on the Carbene Hydroxymethylene”, J. Chem. Educ., 88, 1667-1671, August 2011.
[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.
[9]   | Jensen F, “Introduction to Computational Chemistry”, 3rd Edition, Wiley, 2017.
[10]   | Sathyamurthy, N. and Joseph, T, “Potential energy surface and molecular reaction dynamics”, J. Chem. Educ., 61, 968-971, November 1984.
[11]   | Carey, F. A., and Sundberg, R .J, “Advanced Organic Chemistry, Part A: Structure and Mechanisms”, Springer; 4th edition, 2006.
[12]   | Laidler, K. J, “Chemical Kinetics”, Pearson; 3rd edition, 1997.
[13]   | Truhlar, D (Ed.), “Potential Energy Surfaces and Dynamics Calculations: for Chemical Reactions and Molecular Energy Transfer”, Springer, 2013.
[14]   | Sheppard, D., Terrell, R. and Henkelman, G, “Optimization methods for finding minimum energy path”, J. Chem. Phys., 128, 134106, April 2008.
[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.
[16]   | Quapp, W, “Chemical Reaction Paths and Calculus of Variations”, Theo. Chem. Acc., 121, 227–237, July 2008.
[17]   | Moss, S. J. and Coady, C. J, “Potential energy surfaces and transition state theory”, J. Chem. Educ., 60, 455, June 1983.
[18]   | Feynman, R, “Feynman lectures on Physics”, vol. 2, sec.19, Pearson, 2012.