Supporting Information

Mathematica Tutorial on Minimum Energy Path Calculation on Müller-Brown Potential Energy Surface

Krishnamohan G. P.1* , Omar H.2, Sreeja T. D.3, Roy K. B.4

[1] Department of Science and Humanities, Mar Baselios College of Engineering and Technology (autonomous), Nalanchira, Trivandrum, Kerala 695015, India.
[2] Department of Chemistry, TKM Arts and Science College, Kollam 691005, India.
[3] Department of Chemistry, Sree Krishna College, Guruvayur, Kerala 680102, India.
[4] Department of Chemistry, Sree Neelakanta Government Sanskrit College, Pattambi, Palakkad, Kerala 679306, India.
*Corresponding author: gpkmohan@yahoo.com

In this Mathematica (version 10) notebook, Müller-Brown potential energy surface (MB-PES) is used to find the (1) stationary points, (2) transition structures (TSs) and (3) steepest descent path which connect between TSs and stationary points. The locations of these points were found by using analytical derivatives. The minimum energy path (MEP) path is calculated by applying a steepest descent (SD) algorithm. A one dimensional reaction profile along the MEP has been calculated and plotted. The Appendix is given to illustrate the properties of the Hessian matrix and for a short note on the mass weighted coordinate system along with a brief note on the use MEP in chemical kinetics.

MB-PES and its parameters

Here we define MB-PES and its parameters. This PES is a well known test system for various numerical optimization algorithms. The function is stored in MBpes[x,y].

A0 = {-200.00, -100.00, -170.00, 15.00}; A = {-1.00, -1.00, -6.50, 0.70}; B = {0.00, 0.00, 11.00, 0.60}; 
C0 = {-10.00, -10.00, -6.50, 0.70};x0 = {1.00, 0.00, -0.50, -1.00}; y0 = {0.00, 0.50, 1.50, 1.00};

MBpes[x_, y_]:=Sum[ A0[[term]] Exp[ A[[term]] (x - x0[[term]])^2 + 
B[[term]] (x - x0[[term]]) (y - y0[[term]]) + C0[[term]] (y - y0[[term]])^2 ], {term, 1, 4} ]

To visualize MB-PES, a contour plot of the function is given below:Note that Plot3D is much faster than ContourPlot3D.

MBsurf=ContourPlot[MBpes[x,y],{x,-1.5,1},{y,-0.2,2},Axes->False,Frame->True,ContourShading->False,Contours->70]

Plot3D[MBpes[x,y],{x, -1.5, 1}, {y,-0.2, 2},MeshFunctions->{#3&},Mesh->50,AxesLabel->{x,y}, ColorFunction->"Rainbow"]
(* ContourPlot3D[z==MBpes[x,y],{x,-1.5,1},{y,-0.2,2},{z,-200,111},ColorFunction->"Pastel",MeshFunctions->{#3&},Mesh->20, AxesLabel->{x,y}] *)

The closed (i.e. analytical) form of the MB-potential can be obtained by:

MBpes[x,y]

Transition structures and minima on the 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 first contour diagram for the details. Similarly, one can easily locate two TSs  and one intermediate structure on the PES and are 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 two simultaneous equations, viz., ((∂MBpes[x,y])/(∂ x)) = ((∂MBpes[x,y])/(∂ y)) = 0. 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 latter use. The exact position of TSs, initial, intermediate and final states are given below:

(1) Positon of First Transition structure (TS1):

dfx[x_,y_]=D[MBpes[x,y],x];     (* Analytical derivative: ((∂ MBpes[x,y])/(∂x)) *)
dfy[x_,y_]=D[MBpes[x,y],y];     (* Analytical derivative: ((∂ MBpes[x,y])/(∂y)) *)
FindRoot[{ dfx[x,y]==0,dfy[x,y]==0},{{x,-0.8},{y,0.6}}] (*  The soultion, where the above derivatives of MBpes[x,y] = 0. *)

    xt=-0.8220015587327321; yt = 0.6243128028148713; (* Solution is saved in xt, yt ...*)

(2) Positon of Second Transition structure (TS2):

FindRoot[{ dfx[x,y]==0,dfy[x,y]==0},{{x,0.21},{y,0.21}}] 
(* Solution is saved in xbt, ybt ...*)   
xbt = 0.212486582000662; ybt = 0.292988325107367;

(3) Positon of Initial structure (Initial):

FindRoot[{ dfx[x,y]==0,dfy[x,y]==0},{{x,-0.6},{y,1.4}}] 
(* Solution is saved in xi, yi ...*)   
xi=-0.5582236346330243; 
yi = 1.4417258418046686;

(4) Positon of Intermediate structure (Intermediate):

FindRoot[{ dfx[x,y]==0,dfy[x,y]==0},{{x,0.0},{y,0.5}}] 
(* Solution is saved in xf, yf ...*)   
xf=-0.05001082299820604; yf = 0.4666941048719721;

(5) Positon of Product structure (Final):

FindRoot[{ dfx[x,y]==0,dfy[x,y]==0},{{x,0.6},{y,0.0}}] 
xbf=0.6234994049308766; ybf = 0.028037758528685654;

All of these points can be now represented over the MB-PES contour diagram. Here, we can choose Blue, Cyan, Red and Green dots for the initial, intermediate, TSs and final points of the reaction, respectively (see below).

ContourPlot[MBpes[x,y],{x,-1.5,1},{y,-0.2,2},Axes->False,Frame->True,ContourShading->False,Contours->40,
Epilog->{Blue,Disk[{-0.558,1.44},0.04] , 
Red,Disk[{-0.822,0.624},0.03],Cyan,Disk[{-0.05,0.47},.03],
Red,Disk[{ 0.213,0.293},0.03],Green,Disk[{0.623,0.028 },.04]} ]

Calculation of the Transition Vector

To find the direction of the descent, we need to find a starting vector, usually known as the transition vector or TV (also known as the imaginary mode, since it corresponds to a negative eigenvalue (See (2) for more details). This vector is obtained by diagonalizing the numerical Hessian matrix (See Appendix B) at the TS. For example, at TS1, the numerical Hessian matrix of the TSs and its eigenvalues and eigenvectors are given by:

Clear[x,y, i, j];
DynamicalMx = D[MBpes[x,y],{{x ,y},2}];
DynamicalMx  /.{{x->xt,y->yt}};
MatrixForm[%] 
Eigensystem[% ]

The above 2×2 matrix gives a negative eigenvalue (or strictly, an imaginary frequency proportional to square root of -750.863) and that gives an imaginary mode of [-0.761396 (Overscript[x, ⋀]) + 0.648287 (Overscript[y, ⋀])]. Note that -1*[-0.761396 (Overscript[x, ⋀]) +0.648287 (Overscript[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 (Overscript[x, ⋀]) + 0.865849 (Overscript[y, ⋀])].

Clear[x,y, i, j];
DynamicalMx = D[MBpes[x,y],{{x ,y},2}];
DynamicalMx  /.{{x->xbt,y->ybt}};
MatrixForm[%] 
Eigensystem[% ]

The below figure shows these TVs and indicates the direction of the descent from the corresponding TSs.

CounterMB:=ContourPlot[MBpes[x,y],{x,-1.5,1},{y,-0.2,2},Axes->False,Frame->True,ContourShading->False,Contours->40,Epilog->{Blue,Disk[{-0.558,1.44},.02] , 
Red,Disk[{-0.822,0.624},.02],Cyan,Disk[{-0.05,0.47},.02],Red,Disk[{ 0.213,0.293},0.02],
Green,Disk[{0.623,0.028 },.02]}];

Aarow1={Arrowheads[Medium],Arrow[{{xt,yt},{xt-0.761*0.11,yt+0.648*0.11}}]} ;(* *0.11 is for scaling*)
Aarow2={Arrowheads[Medium],Arrow[{{xt,yt},{xt+0.761*0.11,yt-0.648*0.11}}]};(* *0.11 is for scaling*)
Barow1={Arrowheads[Medium],Arrow[{{xbt,ybt},{xbt-0.500*0.11,ybt+0.866*0.11}}]} ;(* *0.11 is for scaling*)
Barow2={Arrowheads[Medium],Arrow[{{xbt,ybt},{xbt+0.500*0.11,ybt-0.866*0.11}}]};(* *0.11 is for scaling*)
TSvector=Graphics[{Red,Aarow1,Aarow2,Barow1,Barow2}];

Show[CounterMB,TSvector]

The Steepest Descent Path

The general equation for the steepest descent algorithm is given as,

(Subscript[R, i+1])=(Subscript[R, i]) -∇((Subscript[R, i])) · h.

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

The next figure shows a MEP towards the initial/reactant side from the TS1:

step=0.0001;
Backwardwalks = 1; 

Pxb[1] = -0.7613963557703369*step+xt; (* moving away from the x position of the Saddle *)
Pyb[1] = +0.6482866568267857*step+yt; (* moving away from the y position of the Saddle *)

For[i=2,i<1000,
{  Pxb[i]=Pxb[i-1]-dfx[Pxb[i-1],Pyb[i-1]]*step; (* See the SD equation above *)
  Pyb[i]=Pyb[i-1]-dfy[Pxb[i-1],Pyb[i-1]]*step;
  i++;
  If[  Abs[dfx[Pxb[i-1],Pyb[i-1]]] ≤ 0.00001 &&  Abs[dfy[Pxb[i-1],Pyb[i-1]]] ≤ 0.00001,  Break[] ]; (* a STOP condition*)
}
]
Backwardwalks=i

mepBackward:=Table[{Pxb[walks],Pyb[walks]},{walks,Backwardwalks-1}] (* saving the SD data*)
MEPtoLeft = Show[CounterMB,Graphics[ {Dashed,Line[mepBackward]}]]

The above procedure can be used to find the MEP path towards the intermediate structure. In this case TV is reversed from that of the previous case. The result is then plotted as shown below.

Forwardwalks = 1; 

Pxf[1] = +0.7613963557703369*step+xt; (* here, -1.*TS is used *)
Pyf[1] = -0.6482866568267857*step+yt; (* here, -1.*TS is used *)


For[i=2,i<1000,
{ Pxf[i]=Pxf[i-1]-dfx[Pxf[i-1],Pyf[i-1]]*step;
 Pyf[i]=Pyf[i-1]-dfy[Pxf[i-1],Pyf[i-1]]*step;
 i++;
 If[  Abs[dfx[Pxf[i-1],Pyf[i-1]]] ≤ 0.00001 &&  Abs[dfy[Pxf[i-1],Pyf[i-1]]] ≤ 0.00001,  Break[] ];
}
]
Forwardwalks=i

mepForward:=Table[{Pxf[walks],Pyf[walks]},{walks,Forwardwalks-1}] (* saving the SD data*)
MEPtoRight=Show[CounterMB,Graphics[ {Dashed,Line[mepForward]}]]

To view the complete MEP, these plots can be combined and is shown below.

Show[CounterMB,Graphics[ {Dashed,Line[mepBackward]}],Graphics[ {Dashed,Line[mepForward]}]]

Now, the above procedure can be easily repeated from the TS2 using its corresponding TV. For example, from TS2 one can get the MEP towards the intermediate by the following calculation.

step=0.0001;
Backwardwalks2 = 1; 

Pxb2[1] = -0.5003062397007775*step+xbt; (* moving away from the x position of the Saddle-2 *)
Pyb2[1] = +0.8658485240020152*step+ybt; (* moving away from the y position of the Saddle-2 *)

For[i=2,i<1000,
{  Pxb2[i]=Pxb2[i-1]-dfx[Pxb2[i-1],Pyb2[i-1]]*step; (* See the SD equation above *)
  Pyb2[i]=Pyb2[i-1]-dfy[Pxb2[i-1],Pyb2[i-1]]*step;
  i++;
  If[  Abs[dfx[Pxb2[i-1],Pyb2[i-1]]] ≤ 0.00001 &&  Abs[dfy[Pxb2[i-1],Pyb2[i-1]]] ≤ 0.00001,  Break[] ]; (* a STOP condition*)
}
]
Backwardwalks2=i
mepBackward2:=Table[{Pxb2[walks],Pyb2[walks]},{walks,Backwardwalks2-1}]

Similarly, the above procedure is reverted by using a reversed TV of the one used in the above procedure. From this we can calculate MEP from TS2 to the final state or the product.

A complete MEP is shown in the second figure.

Forwardwalks2 = 1; 

Pxf2[1] = +0.5003062397007775*step+xbt; (* here, -1.*TS is used *)
Pyf2[1] = -0.8658485240020152*step+ybt; (* here, -1.*TS is used *)


For[i=2,i<1000,
{ Pxf2[i]=Pxf2[i-1]-dfx[Pxf2[i-1],Pyf2[i-1]]*step;
 Pyf2[i]=Pyf2[i-1]-dfy[Pxf2[i-1],Pyf2[i-1]]*step;
 i++;
 If[  Abs[dfx[Pxf2[i-1],Pyf2[i-1]]] ≤ 0.00001 &&  Abs[dfy[Pxf2[i-1],Pyf2[i-1]]] ≤ 0.00001,  Break[] ];
}
]
Forwardwalks2=i

mepForward2:=Table[{Pxf2[walks],Pyf2[walks]},{walks,Forwardwalks2-1}] (* saving the SD data*)
MEPtoRight2=Show[CounterMB,Graphics[ {Dashed,Line[mepBackward2]}],Graphics[ {Dashed,Line[mepForward2]}]]
Show[%,Graphics[ {Dashed,Line[mepBackward]}],Graphics[ {Dashed,Line[mepForward]}]]

The Reaction Profile

Since we have a full SD path and energies, a 1-D reaction profile can be plotted. Here we chose a reaction coordinate that originates at the TS and is given by,

(Subscript[r, i])=(Underoverscript[∑ , ( i=1), n])(Subscript[r, i-1]) -(Sqrt[(Subscript[Δr, i])]) ,

where (Subscript[r, i]) is the [x, y] coordinate of the MEP path after the i'th SD step, and Δ(Subscript[r, i]) is the Cartesian distance vector (Euclidian norm) from the two consecutive points on the path.

By chosing TS1 as a reference point, we can divide the reaction profile into two parts; (1) towards the reactant and (2) towards the product which passes through the intermediate, TS2.

The following code calculates the 1-D reaction profile towards the reactant from the TS1.

Pxb[0]=xt;Pyb[0]=yt;PES[0]=MBpes[xt,yt];RCb[0]=0;
For[i=1,i<Backwardwalks,{RCb[i]=RCb[i-1]-Sqrt[(Pxb[i]-Pxb[i-1])^2+(Pyb[i]-Pyb[i-1])^2],
PES[i]=MBpes[Pxb[i],Pyb[i]],i++}]
PEStoL = ListPlot[Table[{RCb[n],PES[n]},{n,Backwardwalks-1}]]

MEP towards the product is given below.

RCdat=Join[mepForward,Reverse[mepBackward2],mepForward2];
Totwalks=Forwardwalks+Backwardwalks2+Forwardwalks2-3

(* note: RCdat[[14,1]] ==> RCdat(x) at 14'th point;RCdat[[14,2]] ==> RCdat(y) at 14'th point;*)
(*ListPlot[RCdat,Joined->True]*)

RCf[1]=Sqrt[(xt-Pxf[1])^2+(yt-Pyf[1])^2];
For[i=2,i<Totwalks-1,{RCf[i]=RCf[i-1]+Sqrt[(RCdat[[i,1]]-RCdat[[i+1,1]])^2+(RCdat[[i,2]]-RCdat[[i+1,2]])^2],
PES[i-1]=MBpes[RCdat[[i-1,1]],RCdat[[i-1,2]]],i++}]

PEStoR = ListPlot[Table[{RCf[n],PES[n]},{n,Totwalks}],PlotRange->All]

Show[PEStoL,PEStoR, PlotRange->Automatic,AxesStyle->Thick,AxesOrigin->{-1.25,-160},AxesLabel->{"Reaction coordinate","Energy"}]

Finally, a complete reaction profile can be generated:

Show[PEStoL,PEStoR, PlotRange->Automatic,AxesStyle->Thick,AxesOrigin->{-1.25,-160},AxesLabel->{"Reaction coordinate","Energy"}]

Appendix A: Mass-weighted coordinate system

More realistic (or ab-initio) IRC calculations one should use the mass weighted coordinates . In this coordinate system, the position and force are transformed as, (Subscript[X, mw]) = (Sqrt[m])x and  (Subscript[F(x), mw]) = ((F(x))/(Sqrt[m])) , where m is the mass of the particle, x and F(x) are the Cartesian coordinate position and the force acting on that direction, respectively.  The mass weighted Hessian is calculated from,

(Subscript[H, mw]) = (({{({{(1/(Sqrt[(Subscript[m, 1])])), 0}, {0, (1/(Sqrt[(Subscript[m, i])]))}}), 0}, {0, (1/(Sqrt[(Subscript[m, 3n])]))}})) (Subscript[H,  3n×3n ]) (({{({{(1/(Sqrt[(Subscript[m, 1])])), 0}, {0, (1/(Sqrt[(Subscript[m, i])]))}}), 0}, {0, (1/(Sqrt[(Subscript[m, 3n])]))}})).

Appendix B: Hessian matrix and its diagonalization

The Hessian ( force constant or dynamical matrix) is basically a second derivative matrix which gives the idea of the curvature of the given surface function. As an example, the Cartesian Hessian matrix (as a function of its energy, E((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[y, 1]),(Subscript[y, 2]),(Subscript[z, 1]),(Subscript[z, 2])), where, (Subscript[r, n]) corresponds to the Cartesian position of the atom, n) for a molecule that consits of 2 atoms is defined as,

(Subscript[H,  6 × 6]) = (({{((((∂ )^2)E)/(∂(Subscript[x, 1])∂(Subscript[x, 1]))), ((((∂ )^2)E)/(∂(Subscript[x, 1])∂(Subscript[x, 2]))), ((((∂ )^2)E)/(∂(Subscript[x, 1])∂(Subscript[y, 1]))), ((((∂ )^2)E)/(∂(Subscript[x, 1])∂(Subscript[y, 2]))), ((((∂ )^2)E)/(∂(Subscript[x, 1])∂(Subscript[z, 1]))), ((((∂ )^2)E)/(∂(Subscript[x, 1])∂(Subscript[z, 2])))}, {((((∂ )^2)E)/(∂(Subscript[x, 2])∂(Subscript[x, 1]))), ((((∂ )^2)E)/(∂(Subscript[x, 2])∂(Subscript[x, 2]))), ((((∂ )^2)E)/(∂(Subscript[x, 2])∂(Subscript[y, 1]))), ((((∂ )^2)E)/(∂(Subscript[x, 2])∂(Subscript[y, 2]))), ((((∂ )^2)E)/(∂(Subscript[x, 2])∂(Subscript[z, 1]))), ((((∂ )^2)E)/(∂(Subscript[x, 2])∂(Subscript[z, 2])))}, {((((∂ )^2)E)/(∂(Subscript[y, 1])∂(Subscript[x, 1]))), ((((∂ )^2)E)/(∂(Subscript[y, 1])∂(Subscript[x, 2]))), ((((∂ )^2)E)/(∂(Subscript[y, 1])∂(Subscript[y, 1]))), ((((∂ )^2)E)/(∂(Subscript[y, 1])∂(Subscript[y, 2]))), ((((∂ )^2)E)/(∂(Subscript[y, 1])∂(Subscript[z, 1]))), ((((∂ )^2)E)/(∂(Subscript[y, 1])∂(Subscript[z, 2])))}, {((((∂ )^2)E)/(∂(Subscript[y, 2])∂(Subscript[x, 1]))), ((((∂ )^2)E)/(∂(Subscript[y, 2])∂(Subscript[x, 2]))), ((((∂ )^2)E)/(∂(Subscript[y, 2])∂(Subscript[y, 1]))), ((((∂ )^2)E)/(∂(Subscript[y, 2])∂(Subscript[y, 2]))), ((((∂ )^2)E)/(∂(Subscript[y, 2])∂(Subscript[y, 1]))), ((((∂ )^2)E)/(∂(Subscript[y, 2])∂(Subscript[z, 2])))}, {((((∂ )^2)E)/(∂(Subscript[z, 1])∂(Subscript[x, 1]))), ((((∂ )^2)E)/(∂(Subscript[z, 1])∂(Subscript[x, 2]))), ((((∂ )^2)E)/(∂(Subscript[z, 1])∂(Subscript[y, 1]))), ((((∂ )^2)E)/(∂(Subscript[z, 1])∂(Subscript[y, 2]))), ((((∂ )^2)E)/(∂(Subscript[z, 1])∂(Subscript[z, 1]))), ((((∂ )^2)E)/(∂(Subscript[z, 1])∂(Subscript[z, 2])))}, {((((∂ )^2)E)/(∂(Subscript[z, 2])∂(Subscript[x, 1]))), ((((∂ )^2)E)/(∂(Subscript[z, 2])∂(Subscript[x, 2]))), ((((∂ )^2)E)/(∂(Subscript[z, 2])∂(Subscript[y, 1]))), ((((∂ )^2)E)/(∂(Subscript[z, 2])∂(Subscript[y, 2]))), ((((∂ )^2)E)/(∂(Subscript[z, 2])∂(Subscript[z, 1]))), ((((∂ )^2)E)/(∂(Subscript[z, 2])∂(Subscript[z, 2])))}})).

In this Appendix, using a simple function, f = Sin(x-2.5)Cos(y-2.5), we will demonstrate the evaluation of the Hessian Matrix at the surface saddle point and its diagonalization. The result shows a single negative eigenvalue (i.e. imaginary frequency) exists at the saddle point and the corresponding eigenvector represents the Transition vector (towards the adjacent minimum).

Clear[x,y, i, j];
surfA[x_,y_] := Sin[x-2.5]Cos[y-2.5]
Plot3D[surfA[x,y], {x,-2.5,2.5},{y,-5,5}, AxesLabel->{x,y}, Mesh->None]

Here, we analytically evaluate the Hessian matrix. The resulting matrix is given by a 2×2 matrix, viz. DynamicalMxA.

DynamicalMxA = D[surfA[x,y],{{x ,y},2}] ;
MatrixForm[DynamicalMxA]

We can find the exact location of the TS on the PES by applying the conditions: ((∂ (Subscript[PES, TS]))/(∂ x)) = ((∂ (Subscript[PES, TS]))/(∂ y)) = 0. - which are the two simultaneous equations. The FindRoot function is used to find the roots of these simultaneous equations.

dfxA[x_,y_]=D[surfA[x,y],x];     (* Analytical derivative: ((∂ surfA[x,y])/(∂ x)) *)
dfyA[x_,y_]=D[surfA[x,y],y];     (* Analytical derivative: ((∂ surfA[x,y])/(∂ y)) *)
FindRoot[{ dfxA[x,y]==0,dfyA[x,y]==0},{{x,-0.5},{y,0.6}}] 
(*  The soultion, where surfA[x,y] = 0. *)

xta=-0.6415926535897944; yta = 0.9292036732051044; (* now solution are saved ...*)

Hessian Matrix at TS can be now calculated. The numerical value of this matrix, (({{0, -1}, {-1, 0}})), is used to find diagonalize to find the corresponding eigenvalues and the eigenvectors. The eigenvector (({{-0.70..}, {-0.70..}})), gives the negative eigenvalue (-1).

DynamicalMxA  /.{{x->xta,y->yta}};
MatrixForm[Chop[%] ]
Eigensystem[%]

Finally, a plot is shown to view the direction of the transition vectors. Note that both (({{-0.70..}, {-0.70..}})) and (({{+0.70..}, {+0.70..}})) represent the negative eigenvalue. Arrows originate from (xta, yta). This clearly shows an initial step of the steepest descent path towards the adjacent minima.

ContourPlot[surfA[x, y], {x, -2.5, 2.5}, {y, -5, 5}, AxesLabel -> {x, y}];
arow1={Arrowheads[Large],Arrow[{{xta,yta},{xta-0.707,yta-0.707}}]};
arow2={Arrowheads[Large],Arrow[{{xta,yta},{xta+0.707,yta+0.707}}]};
TSvector=Graphics[{Red,arow1,arow2}];
Show[%%%%,%]

_______________