The reaction path intrinsic reaction coordinate method and the ...

THE JOURNAL OF CHEMICAL PHYSICS 122, 234105 共2005兲

The reaction path intrinsic reaction coordinate method and the Hamilton–Jacobi theory Ramon Crehueta兲 Departament de Química Orgànica Biològica, Institut de Investigacions Químiques i Ambientals de Barcelona, IIQAB-CSIC, Jordi Girona 18, 08034 Barcelona, Catalonia, Spain

Josep Maria Bofillb兲 Departament de Química Orgànica, Universitat de Barcelona, Martí i Franquès 1, 08028 Barcelona, Catalonia, Spain and Centre Especial de Recerca en Química Teòrica, Parc Científic de Barcelona, Baldiri i Reixac 10-12, 08028 Barcelona, Catalonia, Spain

共Received 17 January 2005; accepted 14 April 2005; published online 21 June 2005兲 The definition and location of an intrinsic reaction coordinate path is of crucial importance in many areas of theoretical chemistry. Differential equations used to define the path hitherto are complemented in this study with a variational principle of Fermat type, as Fukui 关Int. J. Quantum Chem., Quantum Chem. Symp. 15, 633 共1981兲兴 reported in a more general form some time ago. This definition is more suitable for problems where initial and final points are given. The variational definition can naturally be recast into a Hamilton–Jacobi equation. The character of the variational solution is studied via the Weierstrass necessary and sufficient conditions. The characterization of the local minima character of the intrinsic reaction coordinate is proved. Such result leads to a numerical algorithm to find intrinsic reaction coordinate paths based on the successive minimizations of the Weierstrass E-function evaluated on a guess curve connecting the initial and final points of the desired path. © 2005 American Institute of Physics. 关DOI: 10.1063/1.1927521兴 I. INTRODUCTION

An interesting problem in theoretical chemistry is to find the lowest energy path of an adiabatic potential energy surface 共PES兲 associated to a molecular rearrangement from a stable configuration state to another, normally associated to reactants and products. This lowest energy path or minimum energy path 共MEP兲 is often referred as reaction path 共RP兲. From a mathematical point of view, the RP is defined as a curved line in the coordinate space, connecting two minima through a first-order saddle point 共FOSP兲, also known as the transition state of a PES. The RPs are the grounds of many dynamical theories. The saddle point gives us the energy barrier, which is the most important quantity to evaluate the transition rates of the rearrangements within the harmonic transition state theory.1 The reformulation of the famous reaction path Hamiltonian 共RPH兲 shows the importance to evaluate first a RP and second to run molecular dynamics constrained to this path. In this way, the RPH methodology and in turn the RP offers the possibility to study dynamics. Using this type of restricted nuclear movement one recovers many important molecular dynamic effects.2 Different definitions of RP produce different curve lines in the PES and the parametrization of these lines, s, is called the reaction coordinate. The reaction coordinate is in fact an arc-length of the RP. In theoretical chemistry, one of the most used RP is the corresponding steepest descent 共SD兲 curve from the FOSP to reactant or product. This SD reaction path in mass weighted coordinates is normally called intrinsic rea兲

Electronic mail: [email protected] Corresponding author: [email protected]

b兲

0021-9606/2005/122共23兲/234105/16/$22.50

action coordinate 共IRC兲.3 However, in some cases the IRC is not the most appropriated MEP.1 The tangent vector of the SD lines and specifically the IRC curve line is characterized by an autonomous system of first-order differential equations. The tangent vector at each point of the line, t关q共s兲兴, is equated to the normalized gradient vector, g共q兲 = ⵱qV共q兲, of this point of the PES, V共q兲,

冏冏

t关q共s兲兴 =

dq共s⬘兲 g关q共s兲兴 , = ds⬘ s⬘=s 兩g关q共s兲兴兩

共1兲

where q共s兲 represents a point of the curve and 兩g关q共s兲兴兩 the norm of the vector g共q兲 evaluated at the point q共s兲 of the curve. The dimension of q共s兲 position vector and the gradient vector g共q兲 is N, N being the number of independent variables. There are a variety of numerical methods to solve the differential equation 共1兲. We propose an analysis of the SD curves and specifically of IRC curve through the calculus of variations and Hamilton–Jacobi theory.4,5 The Hamilton– Jacobi theory has been used in other fields to deal with the integration of ordinary differential equations. This theory formulates the calculus of variations through a nonlinear first-order partial differential equation, called Hamilton– Jacobi equation.6 The integration of the partial differential equations is in general much more difficult than the integration of the ordinary differential equations. The Hamilton–Jacobi theory achieves, in a beautiful manner, to reverse this relationship in some situations. Very often, many applied mathematical problems lead to a system of ordinary differential equations as in the calculus of variations. These equations may be difficult to integrate by elementary methods, but the corre-

122, 234105-1

© 2005 American Institute of Physics

234105-2

R. Crehuet and J. M. Bofill

sponding partial differential equation can be transformed easily into a complete integral by the separation of variables. If the complete integral is obtained, then one solves the corresponding system of characteristic differential equations by a process based on differentiation and elimination.6 Because the proposed IRC reaction path is described mathematically by the total differential equation 共1兲, it should be interesting to investigate the possibility to derive the corresponding Hamilton–Jacobi equation for this type of RP. In addition the resulting Hamilton–Jacobi equation may open other ways to evaluate the IRC curve line. Tachibana and Fukui were the first authors to investigate the IRC path nature from the calculus of variations.7–11 These authors obtained the variational results by using different ways, namely, the Synge’s geodesic principle, the Hamilton principle through an extended Lagrangian, and finally geodesic lines in a Riemannian space.11 The results obtained by these authors have been developed to propose new methods to evaluate the IRC path. The first attempt to find an IRC based on variational methods, as far as we know, is due to Stachó and Bán.12 In fact the algorithm proposed by these authors consists in the variation of a guess curve and the successive correction of this curve based on Mezey’s theory of catchment regions of the gradient field of the PES.13 Other methods to obtain the IRC path based on variations or shifts of an initial guess curve have been proposed and are now widely used.14–24 The IRC path defined by Eq. 共1兲 is in fact an orthogonal trajectory to the contour surface, V共q兲 = const. In the determination of this type of paths the relation between the gradient fields and the associated maps of these orthogonal trajectories is relevant.25 Due to this relation, there exists thus both theoretical and practical reasons for expecting to gain insight into the structural characteristics of the gradient fields of potential energy surfaces and the IRC paths. On the other hand, the basic and complete picture of the problems in the calculus of variations from the Hamilton–Jacobi theory consists in a relation between contours of a surface and curves. These curves are never tangent to these contour surfaces. In addition the difference between two contour lines of these surfaces is related to a functional depending on some arguments that characterize these curves.26 The parallelism between the IRC curve model and the Hamilton–Jacobi theory of the calculus of variations opens the question in both the reformulation of the IRC path from this point of view and new possibilities for determining this type of paths. These two points are the aim of this paper. The SD curves emerging from a stationary point character minimum in the PES can be seen as traveling in an orthogonal manner through the contour lines of this PES. In addition, it should be noted that the construction of contour potential energy surfaces, V共q兲 = const, such that all points satisfying this equation possess the same equipotential difference with respect to another contour line and specifically with the value of the PES in the minimum, is similar to the construction of the Fermat–Huyghens of propagation of the cone rays. This similarity gives the possibility to study the SD curves from Hamilton–Jacobi theory and the calculus of variations point of view. It is important to remember that the Hamilton–Jacobi theory was formulated for the first time to

J. Chem. Phys. 122, 234105 共2005兲

describe, from a rigorous mathematical point of view, the Fermat–Huyghens construction of propagation of light.27 As pointed out before, Tachibana and Fukui studied and analyzed the IRC curve using the calculus of variations and Synge’s principle,10,11 however, we revisit this study starting from a functional like the one that appears in the description of Fermat-type variational principles. We emphasize that the main aim of this work is to formulate, in a consistent and mathematically precise way, the variational character of the IRC curves in terms of the Hamilton–Jacobi theory. Thus we give a sound basis to many of the ideas proposed and used in previous works,24 as well as suggest a new algorithm, which we will describe later. The present work is organized as follows: first, in Sec. II we present a derivation of the SD curves based on the Fermat variational principle. The goal of this section is also to provide the necessary background on the Hamilton–Jacobi theory for a sufficient understanding of this paper. The formulation is written using the symbolism normally employed in the Hamilton–Jacobi theory. We derive the corresponding Hamilton–Jacobi equation, the characteristic equations for the SD curve lines. The concept of the field of extremals applied to SD curves is introduced, too. Two different derivations of Hamilton–Jacobi equation are reported. The second-order necessary and sufficient variational conditions are also studied to analyze the character of the extremal curves, the SD curves. This result is done through the socalled Weierstrass necessary and sufficient conditions and the Jacobi equation.26 In particular, the solution of the Jacobi equation for the SD curve situated in a region sufficiently close to a first-order saddle point is reported and analyzed. In Sec. III a new method to evaluate the IRC path is formulated. In order to obtain a variational method, which actually specifies how to correct the set of parameters that characterize a trial curve, one must expand the corresponding functional beyond first order in the set of variational parameters. However, to obtain sufficient conditions one must allow each parameter to vary in the correct direction. As a consequence, the proposed method is based on the nature of the secondorder necessary and sufficient variational conditions associated with the IRC curve line. Finally an example is given and some conclusions are drawn.

II. THE FERMAT VARIATIONAL PRINCIPLE AS A THEORETICAL BASIS OF IRC MODEL A. The first-order necessary condition

We consider the SD as a path that propagates in a media, the PES, employing the minimum travel potential energy. The travel potential energy of this steepest descent is from a fixed point qR to a variable end point q in the PES. For future purpose, the fixed point qR is taken as a stationary point character minimum of the PES. Using the calculus of variations, the above problem with these conditions can be formulated as follows: we are choosing the curve from the set of all smooth curves q共s兲 all starting from the point qR = q共0兲 and passing through the point q = q共s兲, with tangent dq / ds that minimizes the integral,

234105-3

J. Chem. Phys. 122, 234105 共2005兲

Reaction path and Hamilton–Jacobi theory

⌬VR→q =

冕

s

0

ds⬘ , v共s⬘兲

where s is the arc-length of the path joining the points qR and q, and v共s⬘兲 is a “velocity” or the continuous slowness model. The magnitude v共s⬘兲 is taken non-negative, v共s⬘兲艌 0. In order to connect this variational principle with the SD curve lines and specifically with the IRC, we take using Eq. 共1兲, v共s⬘兲 = 兩g关q共s⬘兲兴兩−1, the inverse of the gradient norm. On the other hand, we consider the arc-length s⬘ to be a function of a new parameter t⬘ defined as ds⬘ =

冑冉冊冉冊 dq dt⬘

T

dq dt⬘ , dt⬘

共3兲

where the superscript T means transposed. With these two definitions, the above variational problem given in Eq. 共2兲 can be reformulated in the form

冕冑冕 t

⌬VR→q共q兲 =

关g共q兲兴T关g共q兲兴冑共dq/dt⬘兲T共dq/dt⬘兲dt⬘

0

F共q,dq/dt⬘兲dt⬘ .

共4兲

0

This integral functional 共4兲 is positive homogeneous of degree one with respect to the tangent vector dq / dt⬘ and does not depend explicitly on the choice of the parameter t⬘ that characterizes the curve.6,26 This type of functional was also proposed some time ago by Fukui et al.10,11 and Elber et al.16 In this formulation of the variational problem, the solution consists in finding the curve q共t⬘兲 connecting both points of the PES, qR = q共0兲 and q = q共t兲, which minimize the integral functional 共4兲. The integral functional ⌬VR→q共q兲 given either in Eq. 共2兲 or Eq. 共4兲 is of the same type as that appearing in Fermat variational principles.6 From the calculus of variations, we apply the basic formula for the general variation of the functional to the problem 共4兲, ⌬VR→q共q兲, regarding the point qR as a fixed point and the point q as a variable. In the region of the PES, ⌬VR→q共q兲 is a single-valued function of the coordinates of the point q. The basic formula which gives the first-order variation of the functional ⌬VR→q共q兲 with respect to q and t is

冕冋 t

␦⌬VR→q =

0

⵱ qF −

d 共⵱dq/dt⬘F兲 = 0. dt⬘

共6兲

If the curve q共t⬘兲 is an extremal, the integral in Eq. 共5兲 vanishes, then the condition ␦⌬VR→q共q兲 = 0 takes the form

兩兵共⵱dq/dt⬘F兲T⌬q + 关F − 共⵱dq/dt⬘F兲Tdq/dt⬘兴⌬t其兩t = 0. 共7兲 Applying these general results to the present problem, the Euler–Lagrange equation 共6兲 takes the form

冢

冉冊冉冊冉冊冉冊

dq dt⬘ I− dq dt⬘

⫻

冢

dq T dt⬘ T dq dt⬘

Hg

冑gTg

−

冣

冑冉冊冉冊冑冉冊冉冊冣 d 2q dt⬘2

dq dt⬘

T

冑gTg

dq dt⬘

dq dt⬘

T

dq dt⬘

= 0,

共8兲

where H is the Hessian matrix, H = ⵱qg , and I is the identity matrix. The dq / dt⬘ vector is the tangent vector of the q共t⬘兲 curve. The solution of Eq. 共8兲 is the autonomous differential equation T

t

=

⵱ qF −

共2兲

d 共⵱dq/dt⬘F兲 dt⬘

册

T

dq = g关q共t⬘兲兴, dt⬘

共9兲

where g关q共t⬘兲兴 is the gradient vector evaluated at the point q共t⬘兲 of the curve. The normalization of Eq. 共9兲 leads to Eq. 共1兲. We conclude that expression 共1兲 is the normalized tangent vector of the path that extremalizes the functional ⌬VR→q共q兲 defined in Eq. 共4兲. In other words, the SD curve connecting the points qR and q is an extremal curve of the variational problem defined in expression 共4兲, ⌬VR→q共q兲. This point may be envisaged in the following way: a path P, starting at the point qR = q共0兲, propagates through the PES according to the speed law or continuous slowness model, v共s兲 = 兩g关q共s兲兴兩−1, arrives at the point q = q共s兲, traveling with the least potential energy variation ⌬VR→q共q兲, as defined in Eq. 共4兲, then this path is characterized by the normalized tangent given in Eq. 共1兲. We note that the condition of least potential energy variation, ⌬VR→q共q兲, as defined in Eq. 共4兲, will be proved below.

z共t⬘兲dt⬘

+ 兵共⵱dq/dt⬘F兲T⌬q + 关F − 共⵱dq/dt⬘F兲Tdq/dt⬘兴⌬t其t , 共5兲 where ⵱Tx = 共⳵ / ⳵x1 , … , ⳵ / ⳵xN兲 , z共t⬘兲 = q*共t⬘兲 − q共t⬘兲 , q*共t⬘兲 and q共t⬘兲 being two neighboring curves in this region of the PES, both starting at the point qR. The curve q*共t⬘兲 connects the points qR , t⬘ = 0 and q + ⌬q , t⬘ = t + ⌬t and the curve q共t⬘兲 connects the points qR , t⬘ = 0 and q , t⬘ = t. The tangent dq / dt⬘ is evaluated at the point q共t兲. Then, the condition ␦⌬VR→q共q兲 = 0 implies that the curve q共t⬘兲 must be an extremal, i.e., a solution of Euler–Lagrange equation 共also known as Euler equation or Lagrange equation兲,

B. Derivation of the Hamilton–Jacobi equation and the corresponding characteristic system of equations

The functional F共q , dq / dt⬘兲, given in Eq. 共4兲, is defined on the curves lying in some region of the PES. We take the unique extremal curve that goes through the point qR to the arbitrary point q. The integral 共4兲 evaluated along this extremal curve, namely, the SD curve, joining these two points takes the value J共q兲 = ⌬VR→q共q兲. Using the language of Hamilton–Jacobi theory, this function J共q兲 is called the geodetic distance between the points qR and q. As explained in the Introduction, to see the IRC method through the Hamilton–Jacobi theory, we use the symbols and the defini-

234105-4

J. Chem. Phys. 122, 234105 共2005兲


tions normally employed in this mathematical theory. Applying Eq. 共5兲 to the extremal curve, SD, at the point q共t兲, and taking the first-order variation in both ⌬q and ⌬t, i.e., ⌬q → dq , ⌬t → dt, and the value of the tangent vector of the extremal curve at this point, dq / dt⬘兩t⬘=t = dq / dt, ␦⌬VR→q共q兲 is transformed into a total differential equation, ␦⌬VR→q共q兲 = dJ共q兲, namely,

␦⌬VR→q共q兲 = dJ共q兲 = 共⵱dq/dtF兲Tdq + 关F − 共⵱dq/dtF兲Tdq/dt兴dt

共10兲

evaluated at the point q , t. From this total differential equation, we introduce the definition of the p vector, namely, ⵱dq/dtF =

dq dt

冑gTg

冑冉冊冉冊 T

dq dt

dq dt

= ⵱ qJ = p

共11兲

and F − 共⵱dq/dtF兲Tdq/dt =

⳵J = 0. ⳵t

共12兲

In the evaluation of both Eqs. 共11兲 and 共12兲, the functional form of F given in the expression 共4兲 has been used. The result of Eq. 共12兲 is the reason why the geodetic distance function J and its total differential form, given in Eq. 共10兲, depend only on q. Proceeding as in the normal calculus of variations, from the expression 共10兲 and using the Legendre transformation, one derives the Hamilton–Jacobi equation.6,26 However in this case the Legendre transformation cannot be applied since the function F is homogeneous of degree one with respect to dq / dt argument. This is the reason why both Eq. 共12兲 and the determinant of the matrix

冉冊冉冊冑冑冉冊冉冊冉冊冉冊 T

T F= ⵱dq/dt⵱dq/dt

g g

dq dt

T

dq dt

冢

dq dt I− dq dt

dq T dt T dq dt

冣

共13兲

vanish. In the present case we proceed in the following way, from Eq. 共11兲 we obtain the dq / dt argument dq = ⵱ qJ dt

冑冉冊冉冊 dq dt

T

冑g T g

dq dt

共14兲

.

From the homogeneity relation of F, we have F=

冉冊 dq dt

T

⵱dq/dtF =

冉冊 dq dt

T

⵱qJ = 冑gTg

冑冉冊冉冊 dq dt

T

dq . dt 共15兲

Multiplying Eq. 共14兲 from the left by 共⵱qJ兲T and using Eq. 共15兲 we obtain 共⵱qJ兲T共⵱qJ兲 = 1. g Tg

共16兲

Equation 共16兲 is a partial differential equation in the ⵱qJ taking the place of the Hamilton–Jacobi equation or eiconal

equation in the present variational problem. The connection between expression 共16兲 and the Hamilton–Jacobi equation is explained in detail in Appendix A, where another derivation of Eq. 共16兲 is also given. As far as we know, this is the first time the Hamilton–Jacobi equation is formulated for the SD path, however, as indicated in the Introduction close results were obtained by Tachibana and Fukui.7–11 Equation 共16兲 tells us that as the parameter t evolves, the coordinates q共t兲 evolve and the contour line with constant potential energy J changes, through the coordinates q, and a point of this contour line is linked to a point of the neighborhood contour line. This set of points defines a curve which extremalizes the functional 共4兲. This curve is in the present case the SD line that goes from the qR point to the q point in the PES. Now we can establish some analogies between the propagation of light through a medium having a variable index of refraction and the present problem. The light rays are given as extremal paths of least time, now the SD curves are extremal paths of the PES. The construction of solutions of the eiconal equation 共16兲 as a contour line with constant potential energy is similar to the Fermat–Huyghens principle for the construction of wave fronts. For future purposes, at this point, we deal briefly with the concept of field of extremal curves which satisfy the trivial Hamilton–Jacobi equation 共16兲. The geodetic distance introduced above is defined from a fixed point qR of the extremal curve to a point of a fixed surface ⌸共q , t兲 = const. This concept of geodetic distance arises by considering the initial point qR of an extremal curve as fixed and seeking a final point q on the given surface ⌸共q , t兲 = const in such a way that the geodetic distance J remains stationary under variations of the point q. Thus in formula 共10兲 we have introduced the value zero for the variation of dq and dt of the end point q since dJ共q兲 = 0 and the initial point qR is fixed. This condition is always satisfied if the point q is varied on the surface ⌸共q , t兲 = const, which means that the vanishing of Eq. 共10兲 is a consequence of the fact that the differentiated form of the surface, ⌸共q , t兲 = const, vanishes, d⌸ = 共⵱q⌸兲Tdq +

⳵⌸ dt = 0. ⳵t

共17兲

Comparing both Eqs. 共10兲 and 共17兲 and using Eqs. 共11兲 and 共12兲 we see that p = ⵱qJ = ⵱q⌸ = ⵱dq/dtF, and ⳵J / ⳵t = ⳵⌸ / ⳵t = F − 共⵱dq/dtF兲Tdq / dt = 0. Notice that in the present case the surface depends only on the coordinates q , ⌸共q兲 = const, and from this we infer that ⌸共q兲 corresponds to J共q兲. The last result is the so-called transversality condition6,26 which in the present problem is a relation between the coordinates q of a point of the surface ⌸共q兲 = const, the p vector, and the tangent dq / dt of the extremal curve, a SD line. Finally in the present problem, the field of extremal curves is centered in the point qR, where all the extremal curves, the SD lines, emerge. Expression 共10兲 gives us the derivatives of J, and its total differential form enables the evaluation of the geodetic distance of an extremal curve as a line integral of this total differential form and also the computation of this geodetic distance independent of the curve used. Let us assume an arbitrary piecewise smooth curve C defined in the region of

234105-5

J. Chem. Phys. 122, 234105 共2005兲


the PES connecting the fixed point qR and the variable end point q. This arbitrary curve, not necessarily an extremal curve 共SD curve兲, is defined at each point of this region by the function qC共t兲 and its tangent by dqC / dt. According to the previous discussion, in the present field of extremals, the characterization of the extremal curves at each point of the considered region, q共t兲, embedded in this field, is given by the corresponding tangent vector dq / dt and the p vector defined in Eq. 共11兲. The function J共q兲, as the line integral of the total differential form defined in Eq. 共10兲 is

冕冋冕冋 q,t

J共q兲 =

共⵱qJ兲TdqC +

qR,0 t

=

共⵱qJ兲T

0

⳵J dt⬘ ⳵ t⬘

册

册

dqC ⳵ J + dt⬘ . dt⬘ ⳵ t⬘

共18兲

Equation 共18兲 is the Hilbert’s invariant integral.6,26 Using Eqs. 共9兲, 共11兲, and 共12兲, Eq. 共18兲 takes the simple form

冕冉 ⬘冊 ⬘ 冕冉 ⬘冊 ⬘ 冕冋冑冑 t

J共q兲 =

pT

dqC dt dt

gCT

dqC dt dt

0

t

=

0

gCT共dqC/dt⬘兲

t

=

0

gCTgC 共dqC/dt⬘兲T共dqC/dt⬘兲

册

F共qC,dqC/dt⬘兲 dt⬘ 共19兲

with p being the vector field at the point qC共t⬘兲 and gC = g关qC共t⬘兲兴. According to the definition of geodetic distance J共q兲 given at the beginning of this section and the expression 共19兲, we have the next equality

冕冑冑冉 ⬘ 冊冉 ⬘ 冊冕冉 ⬘冊 ⬘ 冕冕冋冑冑 ⬘ ⬘ t

J共q兲 =

g Tg

0

t

=

gCT

0

t

=

0

dq dt

dqC dt = dt

T

dq dt

T

dt⬘

t

F共q,dq/dt⬘兲dt⬘

0

gCT共dqC/dt 兲

gCTgC 共dqC/dt 兲T共dqC/dt⬘兲

册

F共qC,dqC/dt⬘兲 dt⬘ , 共20兲

where g is evaluated at the point q共t⬘兲 of the extremal curve. We emphsize that dq / dt⬘ is the tangent vector of the extremal curve, namely, the SD curve joining the points qR and q, and the vector dqC / dt⬘ is the tangent vector of an arbitrary curve joining the same points. Both curves are embedded in the field of the vectors p. In the present variational problem, this field of p vectors are the gradient vectors, p = g共q兲. The two field vectors, namely, the tangent vector of the extremal curves dq / dt⬘ and the p vector, are regarded as a given function of the coordinates q. Equation 共20兲 will be used in the following section. Now, the two vectors, dq / dt⬘ and p, evolve in the field of extremal curves through q共t⬘兲 according to the following equations. From Eqs. 共3兲, 共11兲, and 共16兲 we have

1

dq dq p = T . = dq dt⬘ ds 冑p p dt⬘


T

共21兲

Equation 共21兲 gives us the evolution of the tangent vector of the extremal curve. The equation for the evolution of the p vector is derived from the set of Eqs. 共3兲, 共4兲, 共6兲, and 共11兲, 1

dp dp Hg = T . = 冑 dt ds ⬘ dq g g dt⬘


T

共22兲

Equations 共21兲 and 共22兲 are the canonical system of differential equations coming from the Euler–Lagrange differential equation 共6兲 or, what is identical, the system of differential equations that characterizes the extremal curves of the Hamilton–Jacobi partial differential equation 共16兲.6 The proof of this assertion is given in Appendix B. Finally we accept that the above results are quite trivial since one recognizes that Eq. 共21兲 is the tangent of the SD curve. However, we rewrite these equations in this manner to emphasize that the theory of Hamilton–Jacobi can not only be used to describe and analyze any SD curve and its field but also because it gives us the basis of the proposed algorithm to locate IRC paths explained below. C. The second-order variational conditions

So far, we have only been studying the necessary firstorder variational conditions of a curve to be extremal of the functional integral proposed in Eq. 共4兲. It is interesting to analyze the necessary and sufficient second-order variational conditions to ensure that the extremal curves characterized by Eq. 共9兲 minimize the functional 共4兲. The second-order variational conditions are related to the weak Legendre condition,6,26 and merely consist in analyzing the value and the corresponding sign of the determinant of the matrix, T F = Fdq/dt,dq/dt along each point of the extremal ⵱dq/dt⵱dq/dt curve, where dq / dt is the tangent of the extremal curve. If this determinant is non-negative, det兩Fdq/dt,dq/dt兩艌 0, for all points of the curve, then this curve is an extremal satisfying the second-order necessary weak Legendre condition and possesses a minimum character in a neighborhood. On the other hand, if Fdq/dt,dq/dt is strictly positive, det兩Fdq/dt,dq/dt兩 ⬎ 0, then we say that the extremal curve satisfies the secondorder necessary weak strengthened Legendre condition with minimum character in a neighborhood. For the present case, the value of this determinant at each point of the extremal curve, the SD line, is obtained by substituting in Eq. 共13兲 the value of the tangent vector given in Eq. 共9兲,

冉

T 兩⵱dq/dt⵱dq/dt F兩 dq =g = I − dt

冊

ggT = Fg,g . g Tg

共23兲

From Eq. 共23兲 we see that det Fg,g = 0. This result means that the SD curve satisfies the necessary weak nonstrengthened Legendre condition and makes the integral functional 共4兲 a minimum with respect to continuous comparison functions, q*共t , ␶兲, with its continuous first derivatives in a neighborhood, i.e., as the parameter ␶ → 0 the curve q*共t , ␶兲 → q共t兲 and the tangent dq*共t , ␶兲 / dt → dq / dt, the SD curve and its

234105-6

J. Chem. Phys. 122, 234105 共2005兲


tangent, respectively. A variation, z共t , ␶兲 = q*共t , ␶兲 − q共t兲, which satisfies the previous two conditions as ␶ → 0, is called weak variation, and it means that the extremal curve q共t兲 is compared with curves that approximate q共t兲 in both slope as well as position. In the present problem, an extremal curve q共t兲 which minimizes the integral given in Eq. 共4兲 with respect to all weak variations is called a weak minimum and as shown above satisfies the necessary weak Legendre condition. The Legendre condition is still not sufficient to guarantee a minimum. To do so, we need to compare the extremal curve with all possible curves. These types of comparisons between an extremal curve with any other type of curves is the basis of the second-order variation strong condition of an extremal also known as the Weierstrass necessary condition.26 The Weierstrass sufficient condition in the present problem is obtained if one ensures that the region of the PES containing the extremal curve under consideration 共SD line兲 is covered by the type of field of extremal curves

C ⌬关⌬VR→q兴C = ⌬VR→q 共q兲 − J共q兲 =

冕再冋冕再冋

1−

0

gCT共dqC/dt⬘兲 gCTgC 共dqC/dt⬘兲T共dqC/dt⬘兲

冑

冑

t

=

0

gCTgC冑共dqC/dt⬘兲T共dqC/dt⬘兲 − gCT共dqC/dt⬘兲兴dt⬘

0

t

=

冕关冑 t

1−

defined in Sec. II B. In other words, the extremal curve, SD line, can be embedded in the field of extremals, see below. With this consideration, the Weierstrass sufficient condition is formulated as follows:26 we compare the values of the integral 共4兲 evaluated on both the extremal curve, a SD line joining the points qR and q共t兲, and an arbitrary curve C joining the same points in its neighborhood. The value of the integral 共4兲 evaluated on this extremal curve is J共q兲 according to the discussion of the preceding section. We denote by C 共q兲 the value of the integral 共4兲 evaluated over the ⌬VR→q curve C. By evaluating J共q兲 as an integral along the path C, as formulated in Eq. 共20兲, the comparison of these values is now reduced to a comparison of the integrands alone. According to Sec. II B, the extremal curve, the SD line, is embedded in the field of extremals where the tangent of the extremals at each point of the region of the PES containing this extremal curve, q共t⬘兲, and covered by this field is denoted by dq / dt⬘ = g关q共t⬘兲兴. With these considerations, the difference between these values is

gCT共dqC/dt⬘兲

冑gCTgC冑共dqC/dt⬘兲T共dqC/dt⬘兲

册冑册

冎

gCTgC冑共dqC/dt⬘兲T共dqC/dt⬘兲 dt⬘

冎

F共qC,dqC/dt⬘兲 dt⬘ =

where the integration is taken along the curve C, the gradient vector gC = g共qC兲 and the vector dqC / dt⬘ are the tangent vector of the field and the tangent vector of the curve C embedded in this field, respectively, and both tangents are evaluated at the point qC, a point of the region of the PES covered by this field. To derive Eq. 共24兲, Eqs. 共4兲, 共20兲, and 共3兲 have been used. Equation 共24兲 was first derived by Tachibana and Fukui.9 The integrand of expression 共24兲,

冕关冑 s

gCTgC − gCT共dqC/ds⬘兲兴ds⬘ ,

共24兲

0

C result means that ⌬VR→q 共q兲艌 J共q兲, and consequently the extremal curve, the SD line, connecting the points qR and q共t兲, is actually a strong minimum. In particular, if we take all admissible curves C different from the SD curve, which means that dqC / dt⬘ ⫽ gC, then we have a proper strong miniC 共q兲 ⬎ J共q兲. mum since in this case ⌬VR→q The Legendre conditions are related with the Weierstrass E-function through the expression

E共qC,gC,dqC/dt⬘兲

= 冑gCTgC冑共dqC/dt⬘兲T共dqC/dt⬘兲 − gCT共dqC/dt⬘兲

冋

= 1−

gCT共dqC/dt⬘兲

冑gCTgC冑共dqC/dt⬘兲T共dqC/dt⬘兲

册

E共qC,gC,dqC/dt⬘兲

F共qC,dqC/dt⬘兲,

T F兴dq/dt=g* 共dqC/dt⬘ − gC兲, = 21 共dqC/dt⬘ − gC兲T关⵱dq/dt ⵱dq/dt C

共26兲

共25兲 which in this case is a function of the arguments qC , gC and dqC / dt⬘, is known as Weierstrass E-function. Since for each point qC of the region of the PES covered by the field of extremals with tangent gC and for all possible values of the tangent vector dqC / dt⬘, the Weierstrass E-function is nonnegative, E共qC , gC , dqC / dt⬘兲艌 0, then from the expression 共24兲 we get ⌬关⌬VR→p兴C 艌 0 for all admissible curves C defined in the region of the PES covered by this field. This

T F is that given in Eq. 共13兲 and where the matrix ⵱dq/dt⵱dq/dt * gC = gC + ␪共dqC / dt⬘ − gC兲, with 0 ⬍ ␪ ⬍ 1. The proof of Eq. 共26兲 is given in Appendix C. Since the matrix T F兴dq/dt=g*C is a projector with respect to the vec关⵱dq/dt⵱dq/dt * tor gC, from Eq. 共13兲 and 共26兲 we see that, E共qC , gC , dqC / dt⬘兲 ⬎ 0, if the vector 共dqC / dt⬘ − gC兲 is orthogonal to the vector gC* , otherwise E共qC , gC , dqC / dt⬘兲 = 0. Finally, using Eq. 共24兲–共26兲, we obtain

234105-7

J. Chem. Phys. 122, 234105 共2005兲


FIG. 1. The main scheme of the variational problem given in Eq. 共4兲. See text for a detailed discussion.

C ⌬关⌬VR→q兴C = ⌬VR→q 共q兲 − J共q兲

冕冕兵 t

=

E共qC,gC,dqC/dt⬘兲dt⬘

0

=

1 2

t

T 共dqC/dt⬘ − gC兲T关⵱dq/dt⵱dq/dt F兴dq/dt=g*

0

⫻共dqC/dt⬘ − gC兲其dt⬘ .

C

共27兲

Then the Weierstrass sufficient condition for an extremal curve to be minimal can be formulated, as Eq. 共27兲, to be positive with the possibility to embed the extremal curve joining the two points, q共0兲 and q共t兲, in the field of extremals. Equation 共27兲 will be used in the following section. The basic point in the theory of sufficient conditions just exposed and used above is the possibility of embedding the extremal curve under consideration in a field. If the endpoints of the extremal curve are not too far, i.e., small t parameter, then it can always be embedded in a field. We remember that a field is defined by the set of extremal curves cutting the surface ⌸ transversally.6 The set of extremal curves emerging from a central point will constitute a field up to its envelope or conjugate points to the central point. The first point at which neighboring extremal curves all starting at the same central point intersect is called conjugate point with respect to the central point. In the envelope or conjugate points the set of extremals does not cut ⌸ transversally. In the present problem, the SD curves, emerging from the point qR, which is a stationary point character minimum, intersect for the first time at the stationary points of the PES character saddle points and maxima. These types of stationary point are the conjugate points qCP with respect to central point qR. The above discussion about second-order sufficient conditions is now reduced to the consideration that in the region of the PES containing the SD curve, the extremal curve q共t⬘兲, emerging from the point qR, is a mini-

mum if in the region 0 ⬍ t⬘ ⬍ t, where qR = q共0兲 and q共t兲 = q, a conjugate point of qR does not exist, and Eq. 共27兲 is positive definite in this interval of integration. If the interval of integration 0 ⬍ t⬘ ⬍ t contains tCP , 0 ⬍ tCP ⬍ t, where q共tCP兲 = qCP is a conjugate point with respect to the central point qR, then the extremal curve q共t⬘兲 is not a minimum provided that the Weierstrass E-function defined in Eq. 共26兲 is positive definite, E共qC , gC , dqC / dt⬘兲艌 0, along the extremal curve.26 The conjugate point qCP can be a stationary point character saddle point or maximum in the PES. However, for saddle points with one negative eigenvalue, known as first-order saddle points 共FOSPs兲, only one SD curve emerging form the central point qR arrives at this type of stationary points. As a consequence, the first-order saddle points are not conjugate points with respect to the central point qR, a stationary point with character minimum in the PES. This result is proved from a rigorous mathematical point of view in Appendix D using the Jacobi equation associated to the variational problem under consideration.26 Finally in Fig. 1 we show the basic scheme of all the concepts just exposed for the present variational problem. The explanation of this figure is the following. A SD curve starting at the point q0 of the PES transverses the contour lines V共q兲 − V共q0兲 = const= con1 and V共q + dq兲 − V共q0兲 = const = con2 at the points q and q + dq, respectively. The normal vector of the transversal line, ⌸, at the point q is p1 = ⵱qJ共q兲 = ⵱qV共q兲 = g1, whereas the normal vector of the transversal line, ⌸, at the point q + dq is p2 = ⵱qJ共q + dq兲 = ⵱qV共q + dq兲 = g2, J being the solution of the eiconal equation 共16兲. These transversal lines are tangent to the contour lines. The difference between the contour lines con1 and con2 is denoted by dJ, which is the infinitesimal geodetic distance evaluated on the SD curve using Eq. 共10兲, between the points q and q + dq. The tangent vector of the SD curve at the point q is dq / dt, and dqC / dt denotes the tangent vector at the point q of another arbitrary curve C passing through this point.

234105-8

J. Chem. Phys. 122, 234105 共2005兲


The geodetic distance dJ is taken as a radius of a circle centered at point q, this circle is known as the geodetic circle. In this circle the equality, dJ = F共q , dq / dt兲dt = 关pT1 p1兴1/2关共dq / dt兲T共dq / dt兲兴1/2dt = F共q , dqC / dt兲dt = 关pT1 p1兴1/2 ⫻关共dqC / dt兲T共dqC / dt兲兴1/2dt, is satisfied. The SD path connecting the points q and q + dq is one of the radial curves of this circle. Notice that dJ ⬎ 0, because dt ⬎ 0 and F共q , dq / dt兲 defined in expression 共4兲 is positive definite and due to Eq. 共21兲, the value of the geodetic distance is dJ = pT1 p1dt. Furthermore, comparing the normal vector of the ⌸ line at the point q + dq, we see that it coincides with the radius vector of the geodetic circle at q + dq, consequently both the contour line con2 and the geodetic circle are tangential to each other at that point. The above construction can be extended to all geodetic circles of radius dJ with centers on the contour line con1. As a consequence of these results we say that the contour line con2 = J共q + dq兲 = con1 + dJ with dJ ⬎ 0 , J共q + dq兲 being a solution of the Hamilton–Jacobi equation 共16兲, is an envelope line of the all possible geodetic circles of radius dJ centered on the contour line con1. In fact the above results only show that the geodetic circles are tangential to the contour line con2. Now we need to prove that the line con2 is an envelope line of these circles, in others words, the geodetic circles lie entirely on one side of the contour line con2 apart from the points of tangential contact with it. We will show that this construction gives us a geometrical meaning to both the Weierstrass E-function, Eq. 共25兲, and the Weierstrass necessary condition, Eq. 共24兲. To prove the above question, we consider an arbitrary curve C connecting the points q and q + dqC such that F共q , dqC / dt兲dt = dJ. As a consequence the point q + dqC lies on the geodetic circle of radius dJ, as shown in Fig. 1, and now we need to prove that this point lies on the same side of the contour line con2 as does the point q. We apply the differential form of the Hilbert’s invariant integral, given in Eq. 共19兲, to the arbitrary curve C joining the points q and q + dqC, resulting in dJC = pT1 dqC = pT1 共dqC / dt兲dt = pT1 共dqC / dt兲 ⫻兵关pT1 p1兴1/2关共dqC / dt兲T共dqC / dt兲兴1/2其−1F共q , dqC / dt兲dt = ␳CF共q , dqC / dt兲dt = ␳CdJ. Clearly the values of ␳C are in the domain −1 艋 ␳C 艋 1. Taking into account the values of ␳C and that con2 − con1 = dJ ⬎ 0, the point q + dqC lies on the contour line con1 + dJC = con1 + ␳CdJ, not shown in the figure, such that con2 艌 con1 + ␳CdJ. Consequently the point q + dqC lies on the same side of the contour line con2 as does the point q. Finally, from Eq. 共25兲 we see that the Weierstrass Efunction takes the form dJ − dJC = F共q , dqC / dt兲dt

− ␳CF共q , dqC / dt兲dt = 共1 − ␳C兲dJ = E共q , p1 , dqC / dt兲dt, which is non-negative everywhere as proved in Sec. II C, specifically after Eq. 共24兲. These results show that for the variational problem under consideration, defined in Eq. 共4兲, where the functional F共q , dq / dt兲 is non-negative, the contour line con2 is an envelope line of the geodetic circles centered on the contour line con1.26 The above construction of solutions as envelopes of the present variational problem is exactly the same as that used in the Fermat–Huyghens principle for the construction of wave fronts in the propagation of light,26 as is already explained in the Sec. II B. III. THE WEIERSTRASS SUFFICIENT CONDITION AS A TOOL TO LOCATE THE IRC CURVE A. Background of the algorithm

The above results about the variational nature of the SD curve open the possibility to use a variety of algorithms to integrate this type of curves and specifically the IRC curve line. We present in this section a way to deal with the calculation of the IRC line. The IRC curve is a SD curve line in mass weighted coordinates connecting two minima in the PES, namely, qR and q P, through a first-order saddle point qFOSP.11 Taking into account both the definition of IRC and the concept of centered field discussed in Sec. II, the IRC path can be seen as a SD curve composed by two SD lines, each one being an extremal curve of one of the fields of extremals centered or emerging from the minima qR and q P, and both ending in a common point qFOSP. The qFOSP point is not a conjugate point either for the SD line emerging from the qR point or for the SD line emerging from the q P point. These two centered fields are identical, because for each field the corresponding p vector of this field is a function of the corresponding geodetic distance J , p = ⵱qJ, as explained in the preceding section, and this function J should satisfy the same Hamilton–Jacobi equation 共16兲. As a consequence for both the centered fields the vector of these fields is p = g共q兲. The geodetic distance from the central point qR 共q P兲 to the variable end point qFOSP is denoted by JR共qFOSP兲关J P共qFOSP兲兴. We propose to use the Weierstrass sufficient condition, discussed in the preceding section, as a way to obtain the IRC curve. Since the IRC curve is composed by two extremal curves, SD curves, each one being an extremal belonging to one of the centered fields and for both fields the field vector is p = g共q兲, then the Weierstrass sufficient condition for the IRC curve given in Eq. 共27兲 becomes

⌬关⌬VR→P兴C = ⌬关⌬VR→qFOSP兴C − ⌬关⌬V P→qFOSP兴C C = 关⌬VR→q

=

冕

tf

0

FOSP

− JR共qFOSP兲兴 − 关⌬VCP→q

E共qC,gC,dqC/dt⬘兲dt⬘ =

冕关冑 tf

0

FOSP

− J P共qFOSP兲兴 =

冕

t

E共qC,gC,dqC/dt⬘兲dt⬘ −

0

gCTgC冑共dqC/dt⬘兲T共dqC/dt⬘兲 − gCT共dqC/dt⬘兲兴dt⬘ ,

冕

t


tf

共28兲

234105-9

J. Chem. Phys. 122, 234105 共2005兲


where t f is the value of the independent variable t⬘ in the origin of the field centered in q P, in other words, q P = q共t f 兲. From Eq. 共28兲 it is clear that ⌬关⌬VR→P兴C 艌 0. According to the discussion of preceding section, the sufficiency condition for the IRC path is achieved since it is a SD curve, q共t⬘兲, connecting the points qR = q共0兲 and q P = q共t f 兲 of the PES if in the interval 0 ⬍ t⬘ ⬍ t f , a conjugate point with respect to both points qR and q P does not exist. This shows clearly that the IRC is an extremal curve of the functional given in Eq. 共4兲 with character strong minimum in a neighborhood. Minimizing the ⌬关⌬VR→P兴C function given in Eq. 共28兲 with respect to the parameters that characterize a given arbitrary C curve connecting the fixed end points qR and q P of the PES, iteratively, we will find the near SD curve to this C curve. The curve C is assumed to satisfy the differential equation dqC共t⬘兲 / dt⬘兩t⬘=t = f关qC共t兲兴, where 0 艋 t 艋 t f and f is a vector of a field vector. This curve C is represented as a polygonal line or a chain line defined in the region of the PES under consideration and connecting the points qR and q P and the vector f = ⌬qC , ⌬qC being the difference vector between two consecutive vertex points of the chain. The minimization of the ⌬关⌬VR→P兴C function has the effect of transforming the curve C into another curve such that the field vector f at each point of this new curve coincides as much as possible with the field vector g, which is the field vector of the field of extremals of the functional given in Eq. 共4兲. This is the basis of the proposed algorithm to find the SD curve connecting the stationary points character minimum of the PES, qR and q P. The resulting SD curve is the IRC path if conjugate points do not exist in this curve, in others words, if the SD curve does not contain stationary points character saddle point with more that one negative eigenvalue, as shown in Appendix D. The basis of the algebraic process of the proposed algorithm is the following: first we approximate the integral 共28兲 by using a set of n + 1 points of an arbitrary curve C, ⌬关⌬VR→P兴C =

冕冕关冑 tf


0

tf

=

gCTgC冑共dqC/dt⬘兲T共dqC/dt⬘兲

0

− gCT共dqC/dt⬘兲兴dt⬘ n−1

⬇

共冑gCiTgCi冑⌬qCiT⌬qCi − gCiT⌬qCi兲兺 i=1

= 兵⌬关⌬VR→P兴C其approx n−1

+

共qC⬘i − qCi兲T关兩⵱q 兺 i=1

Ci

Eapprox共qCi,gCi,⌬qCi兲兩C兴 ,

共30兲

is other arbitrary curve and where C⬘ ⵱qCiEapprox共qCi , gCi , ⌬qCi兲 is the Weierstrass E-gradient vector evaluated in the position vectors that characterize the curve C. The explicit form of the Weierstrass E-gradient vector is ⵱qCiEapprox共qCi,gCi,⌬qCi兲 T ⌬qCiHCi = 冑⌬qCi

−

冋冑

冉冑冉冑

gCi

T gCi gCi

T gCi−1 gCi−1

T − 冑gCi gCi

冉冑

−

gCi−1 T gCi−1gCi−1

gCi

T gCi gCi

−

⌬qCi

冑⌬qCiT⌬qCi −

冊

⌬qCi−1

T 冑⌬qCi−1 ⌬qCi−1

⌬qCi

冑⌬qCiT⌬qCi

冊册

,

冊共31兲

where the matrix HCi = H共qCi兲 is the Hessian matrix at the point qCi. The first term of the right-hand side of Eq. 共31兲 involves a Hessian matrix, however, this term can be simplified because the curve C is embedded in a gradient vector field, the field of extremal curves corresponding to the variational problem given in the expression 共4兲. Taking into account this consideration, the term HCi⌬qCi is transformed into HCi⌬qCi = HCi共qCi+1 − qCi兲 ⬇ ⌬gCi = gCi+1 − gCi, where the equality is achieved when the point qCi+1 of the curve C is within the quadratic expansion of the PES centered in qCi. The other term, HCigCi, is the gradient variation vector along the SD curve running through the point qCi, see Eq. 共22兲. As a consequence, HCigCi ⬇ 共g*i − gCi兲 / ␶i, g*i being the gradient vector of a point close to the point qCi and situated on the SD curve that links this point and the qCi point and ␶i the distance between these two points. With these two approximations Eq. 共31兲 is transformed into the expression

⬇

Eapprox共qCi,gCi,⌬qCi兲兺 i=1

= 兵⌬关⌬VR→P兴C其approx ,

兵⌬关⌬VR→P兴C⬘其approx

⵱qCiEapprox共qCi,gCi,⌬qCi兲

n−1

=

this curve, respectively, and their gradient vectors are zero, gR = g共qC0兲 = g P = g共qCn兲 = 0. Now, we expand the function 兵⌬关⌬VR→P兴C其approx up to first order in a Taylor series with respect to the parameters that characterize the arbitrary curve C,

共29兲

n denotes a set of n + 1 position vectors of the where 兵qCi其i=0 arbitrary C curve, ⌬qCi = qCi+1 − qCi, and gCi = g共qCi兲 is the gradient vector evaluated a the point qCi. Notice that this set of selected position vectors is the point vertices of the chain that represents the curve C. The point vectors qC0 = qR and qCn = q P do not appear in the evaluation of 兵⌬关⌬VR→P兴C其approx function, because they are the fixed initial and final points of

冑⌬qCiT⌬qCi 共g*i − gCi兲 − ⌬gCi 冑gCiTgCi ␶i

−

冋冑

T gCi−1 gCi−1

T − 冑gCi gCi

冉冑

冉冑

gCi−1 T gCi−1gCi−1

gCi

T gCi gCi

−

−

⌬qCi

⌬qCi−1

T 冑⌬qCi−1 ⌬qCi−1

冑⌬qCiT⌬qCi

冊册

.

冊共32兲

Notice that in Eq. 共32兲 only gradient and position vectors are involved. Either Eq. 共31兲 or Eq. 共32兲 is the basic expressions for the proposed algorithm. Equation 共32兲 should be used for large molecular systems.

234105-10

J. Chem. Phys. 122, 234105 共2005兲


Finally we have to impose the sufficiency conditions on the algorithm. These sufficiency conditions consist in avoiding the existence of higher order saddle points and maxima in the region of the PES where the search to locate the IRC path is focused. Since these types of stationary points possess higher potential energy with respect to first-order saddle points, during the location process the condition of lower energy at each step is imposed. This restriction is imposed by assuring that at each point vertex of the new C⬘ curve the corresponding value of the potential energy is lower than or equal to the corresponding point vertex of the previous C curve. B. Algorithm description

The minimization scheme of the function 兵⌬关⌬VR→P兴C其approx is carried out by using a Newton–Raphson method solved in a Krylov subspace and its complementary subspace. This method has been employed by Brooks et al.28,29 for both finding minima and locating the IRC curve in the PES. However, there is an important difference with respect to the algorithm proposed by these authors, namely, their objective function is different from that used in the present algorithm which is 兵⌬关⌬VR→P兴C其approx defined in Eq. 共29兲. We recall that this function is related through the calculus of variations with the exact definition of the SD path given in Eq. 共1兲. This function and its derivatives with respect to the position coordinates, given in the expressions either 共31兲 or 共32兲, are well defined. In addition, using this function and the corresponding derivatives we avoid complicated iterative numerical procedures as those that appear in the algorithm described by the above authors.29 The grounds and basic equations of the present method are explained in Appendix E. Now we outline the algorithm. At the initial iteration, ␮ = 1, a guess curve C0 is defined by a set of n + 1 point vectors. Normally, this C0 curve is a straight line. The vectors q0 and qn correspond to the stationary points character minima in the PES related with the geometry structure associated with the reactants and products, respectively. These two vectors are fixed during all the optin−1 , corremization processes. The rest of n − 1 points, 兵qC0i其i=1 spond to the vertices of the chain representation of the curve C0 and are selected in such a way that they are equidistant. The Weierstrass E-gradient vector, ⵱qCiEapprox共qCi , gCi , ⌬qCi兲兩C=C0, for each point vertex, i = 1 , … , n − 1, is computed by using either Eq. 共31兲 or Eq. 共32兲. Taking these gradient vectors as direction vectors, a new set of n − 1 point vertex is computed through the equation qC1i = qC0i + ␩i关⵱qCiEapprox共qCi,gCi,⌬qCi兲兩C=C0兴 ∀ i,i = 1,…,n − 1,

共33兲

where the vector qC1i corresponds to the new position vector of the vertex i of the new curve C1, and the scale factor ␩i is selected in such a way that V共qC1i兲 ⬍ V共qC0i兲. After the evaluation of all new vertex points a reparametrization of the new curve C1 may be necessary to ensure that the vertex points do not cluster. A very simple scheme has been adopted here: when the length of two consecutive segments of this new

curve, ⌬qC1i and ⌬qC1i+1, have a ratio larger than a given threshold, say 0.7, the central point qC1i+1 is moved along its tangent to center it. Finally, the Weierstrass E-gradient vector ⵱qCiEapprox共qCi , gCi , ⌬qCi兲兩C=C1 for each vertex point, i = 1 , … , n − 1, of the new curve C1 is computed. If the convergence criteria is not satisfied for each vertex point, then the rectangular matrices, R共1兲 and G共1兲 for i = 1 , … , n − 1, are i i built and stored according to the expressions 共E2兲 and 共E7兲, respectively, and iteration 2 begins, otherwise the converged curve C1 is the chain representation of the IRC path between the points qR and q P. At the ␮th iteration, for ␮ ⬎ 1, the Newton– Raphson method described in Appendix E is applied. In this case and for each vertex point i, we have, the position vector qC␮−1i, the Weierstrass E-gradient vector ⵱qCiEapprox共qCi , gCi , ⌬qCi兲兩C=C␮−1, and the rectangular matrices R共i ␮−1兲 and G共i ␮−1兲. The set of equations 共E11兲, 共E5兲, and 共E12兲 is solved and the new position vector qC␮i of the new curve C␮ is then computed by a slight modification of Eq. 共E13兲, qC␮i = qC␮−1i + ␩i关A共i ␮−1兲共qC␮i − qC␮−1i兲 + B共i ␮−1兲共qC␮i − qC␮−1i兲兴 ,

共34兲

where the factor ␩i plays the same role as explained in the iteration ␮ = 1. After the evaluation of all new vertex points, a reparametrization of the new curve C␮ may be necessary to ensure that vertex points do not cluster and this is done by using the same procedure as reported in the iteration ␮ = 1. The set of Weierstrass E-gradient vectors, ⵱qCiEapprox共qCi , gCi , ⌬qCi兲兩C=C␮, for i = 1 , … , n − 1, of the new curve C␮ is computed. If the convergence criteria is not satisfied for each vertex point then we update and store the new rectangular matrices, R共i ␮兲 and G共i ␮兲 for i = 1 , … , n − 1, and the new iteration 共␮ + 1兲th begins, otherwise the converged curve C␮ is the chain representation of the IRC path between the points qR and q P. To introduce stability in this minimization algorithm just described, in the evaluation of the Weierstrass E-gradient vector, ⵱qCiEapprox共qCi , gCi , ⌬qCi兲兩C=C␮, using either Eq. 共31兲 or Eq. 共32兲, the current gradient vector gC␮i of the point vortex qC␮i is replaced by the vector 共gC␮i+1 + gC␮i兲 / 2, which is the average gradient vector of the segment vector ⌬qC␮i = qC␮i+1 − qC␮i. In addition, if the angle between the resulting averaged gradient vector gC␮i and the vector ⌬qC␮i is outside the range 关␲ / 2 , −␲ / 2兴, then the sign of ⌬qC␮i vector has to be changed accordingly. In the present implementation of the algorithm, the Weierstrass E-Hessian matrix, ⵱qCi⵱qT Eapprox共qCi , gCi , ⌬qCi兲兩C=C␮, is taken as the unit matrix. Ci

C. Locating the IRC curve on a symmetric potential energy surface

The above proposed integration technique has been applied to find the IRC curve of the surface equation, V共x , y兲 = 2共x2 − 2兲2 + 关共x − y兲2 − 1兴2 + 4共1 − x2 − y 2兲2 + 关共x + y兲2 − 1兴2, which is shown in Fig. 2. The gray arrows are the gradient vectors of the field. As explained in the previous sections,

234105-11


J. Chem. Phys. 122, 234105 共2005兲

FIG. 3. Evolution of the algorithm to find the IRC curve line based on the minimization of the Weierstrass E-function, as defined in Eq. 共29兲, for the PES given in Fig. 2. The solid lines are some contours of the potential energy. The white open dots are the set of 21 points of the initial guess curve. The dark dots indicate the final converged position of the 21 points. In this final position, all points are located in the IRC curve. The IRC curve line follows the sequence of stationary points: R → FOSP→ P. The gradient faced points indicate the behavior of the algorithm during the minimization process.

FIG. 2. Representation of the PES of Sec. III C. The solid lines are some contours of the potential energy. The light solid arrows are some selected gradients of the field, which are the vectors of the field p of the proposed variational problem defined in Eq. 共4兲. The dark dots are the set of 21 points of the initial curve. The point R is labeled as 1 and the point P as 21. The bold faced arrows are the Weierstrass E-gradient vectors ⵱qiEapprox共qi , gi , ⌬qi兲兩C共R→P兲 associated with the points of the initial curve.

these gradient vectors are the tangent vectors of the SD curves and these curves define the field of extremals. Due to the symmetry of the surface, only two SD paths are the IRC curves. These two paths follow the sequence of stationary points R → FOSP→ P, and correspond to the SD paths connecting R and P, such that the integral given in Eq. 共4兲 evaluated on this curve has the lowest value. This surface is challenging because it presents a maximum, which is an attractor of the SD lines. The algorithm has to make sure that it converges to the SD line that passes through FOSP. For the sake of simplicity, in this section we have dropped the subscript C in the qi , gi, and ⌬qi vectors. The set of open dots defines a broken line, which is the initial guess curve, C共R → P兲. This initial curve is characterized by n = 21 points, where the points q1 = qR and q21 = q P are fixed and the rest of the points are allowed to move. The bold faced arrows correspond to the set of 19 vectors, ⵱qiEapprox共qi , gi , ⌬qi兲兩C共R→P兲 for i = 2,…, 20, of the initial curve C共R → P兲. The direction of these vectors is different depending on the point, and for points 9 and 13 a decrease of the Weierstrass E-function would imply an increase of potential energy. Therefore the ␩i factors applied to the points where gTi 关⵱qiEapprox共qi , gi , ⌬qi兲兩C共R→P兲兴 ⬍ 0 for i = 2,…, 20 are

chosen positive and vice versa. In this manner the new set of generated points will possess lower potential energy. Finally, the algorithm converges in the way that all 21 points are located in the IRC curve. This final position is represented by a set of dark dots in Fig. 3. As is conventional in many minimization procedures, we take some steps as steepest descent, until we get close to the quadratic region, by using Eq. 共33兲. In this example, 10 steps

were taken this way. This procedure also allows to generate a set of matrices, R共i ␮兲 and G共i ␮兲, as defined in Eqs. 共E2兲 and 共E7兲, respectively, for i = 2,…, 20. At every Newton– Raphson step and for each qi point, we first evaluate the new position due to the space spanned by the set of difference positions using Eqs. 共E11兲 and 共E5兲. If the new point implies a descent V共x , y兲, then it is accepted, otherwise the Newton– Raphson step is rejected and a steepest descent step is taken. Finally the variation due to the complementary subspace is computed through Eq. 共E12兲 and properly scaled by ␩i to guarantee a descent effect in the function Eapprox共qi , gi , ⌬qi兲兩C共R→P兲. For the sake of completeness, we have compared the present method with the nudged elastic band method17,18 共NEB兲 with a tangent vector defined as ⌬qi = qi+1 − qi for i = 2,…, 20 on the current curve C共R → P兲. In such a case the algorithm is unstable and, as expected, leads to kinks.18 We are aware that corrections to this are possible,18 but this comparison was only done to show that the method presented here does not suffer from this problem and is stable even with a crude tangent estimation. The convergence of this method has also been tested. Figure 4 shows the decrease in the function 兵⌬关⌬VR→P兴C其approx defined in Eq. 共29兲, being C共R → P兲, the current curve, for the Newton–Raphson algorithm described here and a quenched velocity Verlet as the one used in NEB.17,18 It has to be mentioned that each quenched velocity Verlet step involves a single gradient evaluation for each point, whereas the Newton–Raphson step may involve several gradient evaluations to choose the correct ␩i of Eq. 共34兲. However, we have seen that the number of gradient evaluations is, on average, close to 1.5 per point, or even less when we get closer to the optimized path. The use of the approximate Weierstrass E-gradient vector, through Eq. 共32兲, is also surprising, because its convergence behavior is excellent until very close to converged IRC path 共Fig. 4兲. At that point, a more accurate gradient would be needed to decrease the value of the function 兵⌬关⌬VR→P兴C其approx. Equation 共共32兲兲 is therefore useful to get

234105-12

J. Chem. Phys. 122, 234105 共2005兲


Hamilton–Jacobi equation associated with this Fermat variational principle has been derived, resulting in a very simple expression. As in the normal calculus of variations, from the derived eiconal equation we obtain the associated characteristic system of equations of the IRC curve. The analysis of the second-order variation permits to establish the strong minimum character of the IRC path and from this analysis to propose new algorithm to locate IRC curves.

ACKNOWLEDGMENTS

FIG. 4. The behavior of log兵⌬关⌬VR→P兴C其approx, defined in Eq. 共29兲, vs iteration number, during the minimization process of the Weierstrass E-function, for the example discussed in Fig. 2 and 3. See the text for more details.

pretty accurate paths that can then be refined by increasing the number of points in the discretization or with conventional transition state searches. It is also worth pointing out that the reparametrization of the current C共R → P兲 is only necessary during the initial steps, because the Weierstrass E-gradient vector, evaluated by using Eq. 共31兲, once we are close to the final SD curve, is perpendicular to this curve. This is also true for the perpendicular force used in the NEB method, but the difference is that Eq. 共31兲 corresponds to a gradient vector of a certain objective function while the NEB projection makes its force nonconservative. Certainly a number of authors have also worked on the convergence of NEB or similar chain-of-states methods such as that used in the proposed algorithm.19–21,23,29 VandenEijden et al.19 and Wales et al.21 implemented a Broyden– Fletcher–Goldfarb–Shanno 共BFGS兲 minimizer type,30 both achieving superlinear convergence. We have also commented briefly on the differences and similarities of the recent NEB improved algorithm described in Ref. 29. The NEB method has also been coupled to Car–Parrinello molecular dynamics23 for better convergence with density functional theory calculations. All of these methods outperform the steepest descent minimization of the original NEB method that we, and the other reported authors, have used for comparison. It would surely be interesting to compare all of them but this is beyond the scope of this work. Indeed, the algorithm presented in this paper was mainly formulated to prove and to show the potential applications of a correct description of IRC method via Hamilton–Jacobi theory. A computationally less expensive formulation of the method would be desirable and is part of our future work, nevertheless, we would like to stress that the described algorithm represents an improvement of the original NEB method. IV. CONCLUSIONS

The calculus of variations and specifically the Fermat variational principle can be used as a tool to study the nature of the IRC model and to establish new algorithms to evaluate this type of path. From this point of view, the IRC paths are extremal curves of a Fermat variational principle. The

Financial support from the Spanish Ministerio de Ciencia y Tecnología 关Project No. CICyT BQU2002-00293 and the Ramón y Cajal program 共R.C.兲兴 and, in part, from the Generalitat de Catalunya 共Project No. 2001SGR-00048兲 is fully acknowledged.

APPENDIX A: ALTERNATIVE DERIVATION OF THE HAMILTON–JACOBI EQUATION FOR THE EXTREMAL CURVES OF THE TYPE SD LINE

In this appendix we establish the connection between Eq. 共16兲 and the Hamilton–Jacobi equation or, in others words, we show that Eq. 共16兲 plays the role of a Hamilton– Jacobi equation. To prove this assertion, the basic idea is to deal with homogeneous functional with respect to the argument dq / dt, but of degree greater than one, because in this case it is possible to apply the Legendre transformation and then to derive in the standard way the corresponding Hamilton–Jacobi equation.6,26 In addition this new homogeneous functional should be related to the functional of expression 共4兲. To this aim, we first choose the function q共x兲 such that variation of the integral I共q兲 =

冕冉冊冕冉冊冉冊 x⬘

S q,

0

dq dx = dx

x⬘

0

g Tg

dq dx

T

dq dx dx

共A1兲

with respect to this function vanishes. The variable x is proportional to the variable t of variational problem 共4兲 in the way that dx = F共q , dq / dt兲dt, where both t and F共q , dq / dt兲 are those given in expression 共4兲. We note that the functional S共q , dq / dx兲 of Eq. 共A1兲 is homogeneous of degree two with respect to the argument dq / dx. The extremal curves q共x兲 of the variational problem 共A1兲 satisfy the Euler–Lagrange differential equation ⵱ qS −

d 共⵱dq/dxS兲 = 0. dx

共A2兲

The tangent of these extremal curves is dq / dx = g / 共gTg兲. Due to the homogeneity of the S functional, it satisfies the relation S=−S+

冉冊 dq dx

T

共⵱dq/dxS兲.

共A3兲

Differentiating Eq. 共A3兲 with respect to x and using Eq. 共A2兲, we conclude that the functional S becomes constant along an extremal curve q共x兲, because dS / dx = 0. In other words, S = const= 1 along an extremal curve. With this result, Eq. 共A2兲 can be transformed into the following way:

234105-13

冑冋 1

2 S

J. Chem. Phys. 122, 234105 共2005兲


⵱ qS − 1

册

1 ⵱dq/dxS dS d ⵱dq/dxS + 4 S dx dx d

冉

1

冊

=

⵱ qS − ⵱dq/dxS dx 2冑S 2 冑S

=

dt d ⵱qF − 共⵱dq/dtF兲 = 0. dx dt

冋

册

APPENDIX B: THE CHARACTERISTIC SYSTEM OF DIFFERENTIAL EQUATIONS OF THE HAMILTON–JACOBI EQUATION „16…

共A4兲

Since, dt / dx = 1 / F共q , dq / dt兲 ⫽ 0, any extremal curve satisfying Eq. 共A2兲 also satisfies Eq. 共6兲. Consequently the Hamilton–Jacobi equation of the variational problem 共A1兲 is related with the Hamilton–Jacobi equation of the variational problem 共4兲. To prove this assertion, first we derive the Hamilton–Jacobi equation for the variational problem 共A1兲. Since S is an homogeneous functional of degree two with respect to the argument dq / dx, or in other words, det T S兴 = 2gTg ⫽ 0, except in stationary points, we 关⵱dq/dx⵱dq/dx can use the Legendre transformation. Defining the vector p* as p* = ⵱dq/dxS = 2gTg

dq dx

共A5兲

and substituting in the right-hand side of Eq. 共A3兲, the expression of the vector dq / dx as a function of p* obtained from Eq. 共A5兲, we get the transformed Legendre function of S, namely, L共q,p*兲 =

1 共p*兲Tp* . 4 g Tg

共A6兲

Let us assume that J共q兲 is the geodetic distance function and a solution of the partial differential equation 共16兲; then we want to characterize the field of extremals such that these extremal curves transverse all families of contour line surfaces, J共q兲 = const, as explained in Sec. II B. Now we define a field vector p in the region of the PES considered by the equation Substituting Eq. 共B1兲 in Eq. 共16兲 we have p Tp = 1. g Tg

共B2兲

Now we use the family of curves defined by the ordinary differential equation 共14兲, 1

dJ*共q,x兲 = 共⵱dq/dxS兲Tdq + 关S − 共⵱dq/dxS兲Tdq/dx兴dx

冉冊

⳵J dx, ⳵x *

= p*dq − L共q,p*兲dx = 共⵱qJ*兲Tdq +

共A7兲 where Eqs. 共A3兲 and 共A5兲 have been used. From Eq. 共A7兲 we have p* = ⵱qJ*共q , x兲 and 1 共⵱qJ*兲T⵱qJ* ⳵ J*共q,x兲 = − L共q,⵱qJ*兲 = − . 4 ⳵x g Tg

共A8兲

In the derivation of expression 共A8兲, expression 共A6兲 has been used. Equation 共A8兲 is the Hamilton–Jacobi formula of the variational problem 共A1兲. Since the right-hand side of Eq. 共A8兲 is independent of x, a solution of this partial differential equation is J*共q , x兲 = Cx + 2J共q兲, where C is a constant. Using the fact that L共q , p*兲 = 1, we have C = −1 and as a result we obtain Eq. 共16兲.

dq dq p , = = dq dt⬘ ds 冑pTp dt⬘


T

共B3兲

where the vectors ⵱qJ共q兲 and g that appear on the right-hand side of Eq. 共14兲 are substituted by the vector p by using Eqs. 共B1兲 and 共B2兲, respectively. Equation 共3兲 has also been used to change the independent variable t⬘. The vector p becomes a function of the independent variable t⬘ along the curve characterized by the system of differential equations 共B3兲. By differentiation of the vector p with respect to this independent variable t⬘, we obtain 1

We emphasize that the left-hand side of Eq. 共A3兲 is equal to the transformed Legendre function L共q , p*兲 and, as a consequence, in this case L共q , p*兲 = S共q , dq / dx兲. If J*共q , x兲 is the geodetic distance corresponding to the variational problem 共A1兲, then its total differential form, whose general formula is given in Eq. 共10兲, now reads

共B1兲

p = ⵱qJ共q兲.

dp dq = 关⵱q⵱qTJ共q兲兴 dt⬘ dq dt⬘ dt⬘


T

1


T

dq dt⬘

.

共B4兲 By differentiation of the partial differential equation 共16兲 with respect to q, and after some trivial rearrangements we obtain the identity

关⵱q⵱qTJ共q兲兴冑

p T

p p

−H

g

冑gTg = 0.

共B5兲

Finally using Eqs. 共B3兲–共B5兲 and 共3兲, we get 1

g dp dp =H T . = 冑g g dq dt⬘ ds dt⬘


T

共B6兲

Thus Eqs. 共B3兲 and 共B6兲 or Eqs. 共21兲 and 共22兲 characterize a family of curves as extremal curves and correspond to the characteristic system of differential equations associated to the partial differential equation 共16兲 and consequently these extremal curves are the extremals of the variational problem given in Eq. 共4兲. APPENDIX C: PROOF OF EQ. „26…

The functional F共q , dq / dt⬘兲 defined in Eq. 共4兲 can be expanded with respect to the argument dq / dt⬘ using a Tay-

234105-14

J. Chem. Phys. 122, 234105 共2005兲


lor’s series up to first order with the remainder around the point dq / dt⬘ = g, the tangent of the extremal curves of the field of extremals, F共qC,dqC/dt⬘兲 = F共qC,gC兲 + 共dqC/dt⬘ − gC兲T⵱dq/dt⬘ ⫻F共qC,dq/dt⬘兲兩dq/dt⬘=gC + 21 共dqC/dt⬘ − g C兲

T

T 关⵱dq/dt⬘⵱dq/dt⬘F共qC,dq/dt⬘兲兴dq/dt⬘=g* C

⫻共dqC/dt⬘ − gC兲,

共C1兲

where the vector gC* = gC + ␪共dqC / dt⬘ − gC兲 with 0 ⬍ ␪ ⬍ 1 , qC is a point of the region of the PES covered by the field of extremal curves, the SD lines emerging from the fixed point qR , dqC / dt⬘ is the tangent vector of an arbitrary curve C embedded in this field of extremals at the point qC and the gradient vector gC = g共qC兲 is the tangent vector of the extremal curve also called tangent of the field at this point qC. Now we rearrange Eq. 共C1兲 as follows: F共qC,dqC/dt⬘兲 − F共qC,gC兲 − 共dqC/dt⬘ − gC兲T⵱dq/dt⬘ ⫻F共qC,dq/dt⬘兲兩dq/dt⬘=gC = 21 共dqC/dt⬘ − gC兲T ⫻关⵱dq/dt⬘⵱dq/dt⬘F共qC,dq/dt⬘兲兴dq/dt⬘=g*共dqC/dt⬘ − gC兲. T

C

共C2兲 Substituting in the left-hand side of Eq. 共C2兲 the expressions for F共qC , dqC / dt⬘兲 and ⵱dq/dt⬘F共qC , dq / dt⬘兲, which are given in the Eqs. 共4兲 and 共11兲, respectively, and after some trivial rearrangements we get

冑gCTgC冑共dqC/dt⬘兲T共dqC/dt⬘兲 − gCT共dqC/dt⬘兲 = 21 共dqC/dt⬘ − gC兲T关⵱dq/dt⬘⵱dq/dt⬘ T

⫻F共qC,dq/dt⬘兲兴dq/dt⬘=g*共dqC/dt⬘ − gC兲. C

共C3兲

Finally, substituting Eq. 共C3兲 in Eq. 共25兲 we obtain Eq. 共26兲. APPENDIX D: ANALYTICAL REPRESENTATION OF CONJUGATE POINTS BELONGING TO SD LINE, PROOF THAT A SD CURVE CONNECTING BOTH A MINIMUM AND A FIRST-ORDER SADDLE POINT IN THE PES AND EMBEDDED IN A CENTERED FIELD OF SD CURVES DOES NOT CONTAIN CONJUGATE POINTS

The origin of the concept of conjugate points of a field of extremal curves is related to the following question: given an extremal curve q共t⬘兲, i.e., a curve satisfying Eq. 共6兲, and the varied curve q*共t⬘兲 = q共t⬘兲 + z共t⬘兲, which conditions have to be imposed on z共t⬘兲 such that the varied curve q*共t⬘兲 is also an extremal curve satisfying Eq. 共6兲? To answer this question first we substitute the curve q*共t⬘兲 = q共t⬘兲 + z共t⬘兲 into the Euler–Lagrange equation 共6兲 ⵱qF共q + z,dq/dt⬘ + dz/dt⬘兲 −

d 关⵱dq/dt⬘F共q + z,dq/dt⬘ + dz/dt⬘兲兴 = 0. dt⬘

共D1兲

Second, taking into account that q共t⬘兲 is also a solution of the Euler–Lagrange equation 共6兲, using Taylor’s series and neglecting infinitesimal order higher than one with respect to

both z共t⬘兲 and dz共t⬘兲 / dt⬘, and finally combining terms, we obtain the linear differential equation

冋

dz共t⬘兲 d T 共⵱dq/dt⬘⵱dq/dt⬘F兲 dt⬘ dt⬘

冋

− 共⵱q⵱qTF兲 −

册

册

d T 共⵱q⵱dq/dt⬘F兲 z共t⬘兲 = 0. dt⬘

共D2兲

Equation 共D2兲 is the Jacobi equation26 of the variational problem given in expression 共4兲. The Jacobi equation, except for infinitesimals of order higher than one with respect to z共t⬘兲 and dz共t⬘兲 / dt⬘, is the linear differential equation satisfied by the difference between two neighboring or infinitely close extremal curves. Given two neighboring extremal curves, q*共t⬘兲 and q共t⬘兲, of a centered field of extremals, emerging from the same initial point, qR = q共0兲, the difference vector function, z共t⬘兲 = q*共t⬘兲 − q共t⬘兲, satisfying Eq. 共D2兲 is a nonzero solution of Jacobi equation within an infinitesimal order higher than one relative to z共t⬘兲 and dz共t⬘兲 / dt⬘. In the present case, N solutions of the Jacobi equation 共D2兲 exist, N being the dimension of q vector. The set of initial conditions for each one of these N solutions of the Jacobi equation 共D2兲 are obtained as follows: the initial point is the central point, q共0兲 = qR, and the vector difference z共0兲 is a set equal to the zeroed vector, z共0兲 = 0, whereas its first derivative with respect to t⬘ , dz共t⬘兲 / dt⬘兩t⬘=0 = 1, where 1T = 共11 , … , 1N兲. Given an extremal curve q共t⬘兲 the point qCP = q共tCP兲 is said to be a conjugate point to the central point of the field qR = q共0兲, if at qCP = q共tCP兲 the difference q*共tCP兲 − q共tCP兲, q*共t⬘兲 being a neighboring extremal curve emerging from the same initial point qR = q共0兲, is an infinitesimal of order higher than one relative to z共tCP兲 and dz共t⬘兲 / dt⬘兩t⬘=tCP. Now, we apply these results to the SD curve emerging from the qR point that arrives to the first-order saddle point qFOSP. Substituting the integrand F共q , dq / dt⬘兲, given in expression 共4兲, into Eq. 共D2兲, after some rearrangement we obtain the corresponding Jacobi equation for the variational problem 共4兲, − gTgP

冉

d2z共t⬘兲 gTHg T 冑共I + Q兲 − HQ − g g dt⬘2 g Tg

冊

− QH

dz共t⬘兲 + 共HPH − HPHQ − HQHP dt⬘

− PHQH − QHPH兲z共t⬘兲 = 0,

共D3兲

where I is the unit matrix, g is the gradient vector, and H is the Hessian matrix at the point q of the SD curve, the matrices P and Q are the projectors, P = I − ggT / 共gTg兲, and Q = ggT / 共gTg兲, respectively. In the derivation of Eq. 共D3兲 only the quadratic terms in the expansion of the PES around the point q have been considered. As a consequence the whole integration of the Jacobi equation from qR = q共0兲 to qFOSP = q共t兲 associated to a SD curve is carried out by stepwise quadratic approximation, using Eq. 共D3兲. Notice that t⬘ = 0 at the qR point and t⬘ = t at the qFOSP point. The vector z共t⬘兲 and its derivatives with respect to t⬘ are expressed as a linear combination of the eigenvectors of the Hessian matrix H,

234105-15

J. Chem. Phys. 122, 234105 共2005兲


namely, z共t⬘兲 = Va共t⬘兲 , dz共t⬘兲 / dt⬘ = Va⬘共t⬘兲, and d2z共t⬘兲 / dt⬘2 = Va⬙共t⬘兲, V being the matrix defined by the set of orthonormal eigenvectors of the H matrix. When the SD curve arrives at the first-order saddle point qFOSP, then around this point in the direction of the SD curve, the gradient vector can be approximated as g ⬇ CvTV, where vTV is the column vector of the V matrix corresponding to the normalized eigenvector of the Hessian matrix with negative eigenvalue and C a small scalar taking the value C = 0 at the first-order saddle point.31 With these considerations, at the first-order saddle point, q共t兲 = qFOSP, the first two terms of Eq. 共D3兲 vanish. On the other hand, since at the first-order saddle point the matrix T T and the matrix Q = vTVvTV , then the terms, P = I − vTVvTV PHQ= QHP= O, O being the zeroed matrix. As a consequence the remainder term of Eq. 共D3兲 is HPHVa共t兲 = 0. Since we are interested in the SD curve that arrives at the first-order saddle point, we multiply from the left the remainT T and I − vTVvTV resulting in 关a共t兲兴TV ing term by both vTVvTV ⫽ 0 and 关a共t兲兴i = 0 ∀ i = 1 , N and i ⫽ TV, respectively. Taking into account the initial conditions, z共0兲 = 0 and dz共t⬘兲 / dt⬘兩t⬘=0 ⫽ 0, and invoking the continuous dependences of the solution of Eq. 共D3兲 with respect to these initial conditions, then a共t兲 ⫽ 0 , a⬘共t兲 ⫽ 0, and a⬙共t兲 ⫽ 0. Consequently, the solution of Eq. 共D3兲 related with the SD curve that arrives at the first-order saddle point, qFOSP = q共t兲, from the is such that z共t兲 minimum point, qR = q共0兲, ⫽ 0 , dz共t⬘兲 / dt⬘兩t⬘=t ⫽ 0 and d2z共t⬘兲 / dt⬘2兩t⬘=t ⫽ 0, by invoking the transformation z共t⬘兲 = Va共t⬘兲 , dz , 共t⬘兲 / dt⬘ = Va⬘共t⬘兲, and d2z共t⬘兲 / dt⬘ = Va⬙共t⬘兲 at t⬘ = t. This result shows that the firstorder saddle point is not a conjugate point with respect to the central point qR, a stationary point character minimum in the PES. APPENDIX E: MATHEMATICAL BASIS OF THE NEWTON–RAPHSON METHOD SOLVED IN BOTH A KRYLOV SUBSPACE AND ITS COMPLEMENTARY SUBSPACE USED IN THE MINIMIZATION OF THE FUNCTION GIVEN IN EQ. „29…

In this appendix we derive the set of equations of the Newton–Raphson algorithm to be applied in the minimization of the function given in Eq. 共29兲 with respect to the parameters that characterize the arbitrary C curve to locate a SD curve. The Newton–Raphson equations are projected and solved in both a Krylov subspace and its complementary subspace. The Krylov subspace is generated during the minimization process. This method has been reviewed several times and used in different contexts.28,29,32–34 At the 共␮ + 1兲th iteration of the minimization process, ␮ ⬎ 1, the current curve is denoted by C␮, the point vector of the vertex i by qC␮i, and the Weierstrass E-gradient vector by ⵱qCiEapprox共qCi , gCi , ⵱ qCi兲兩C=C␮. The dimension of these two vectors is N. The Weierstrass E-gradient vector can be evaluated by either Eq. 共31兲 or Eq. 共32兲. There exists a set of ␮ vectors 兵qCvi其␮v=1 and the corresponding Weierstrass E-gradients 兵⵱qCiEapprox共qCi , gCi , ⌬qCi兲兩C=Cv其␮v=1 for each point vertex i. We define the matrix A共i ␮兲 which is a projector onto the subspace defined by the vector differences 兵共qCvi兲 −1 . The matrix B共i ␮兲 corresponds to the projector onto − qC␮i其␮v=1

the complementary subspace. Both A共i ␮兲 and B共i ␮兲 are the matrices of dimension N ⫻ N. These matrices have the following properties: A共i ␮兲 + B共i ␮兲 = I , A共i ␮兲A共i ␮兲 = A共i ␮兲 , B共i ␮兲B共i ␮兲 = B共i ␮兲, and A共i ␮兲B共i ␮兲 = O, where I is the N ⫻ N identity matrix and O is the N ⫻ N zeroed matrix. The explicit form of the A共i ␮兲 matrix is A共i ␮兲 = R共i ␮兲关共R共i ␮兲兲T共R共i ␮兲兲兴−1共R共i ␮兲兲T ,

共E1兲

R共i ␮兲

where is a rectangular matrix of dimension N ⫻ 共␮ − 1兲 whose column vectors are defined by vector differences between the current and the previous position vectors, R共i ␮兲 = 关共qC1i − qC␮i兲,…,共qC␮−1i − qC␮i兲兴 .

共E2兲

The set of column vectors that defines the R共i ␮兲 rectangular matrix is assumed to be linearly independent. The application of the Newton–Raphson method to minimize the function 兵⌬关⌬VR→P兴C其approx, defined in Eq. 共29兲, at the point i for the iteration ␮, results in the following expression:

关兩⵱

兩

T qCi⵱qCiEapprox共qCi,gCi,⌬qCi兲 C=C ␮

兴共q

C␮+1i

− q C␮i兲

= 兩− ⵱qCiEapprox共qCi,gCi,⌬qCi兲兩C=C ,

共E3兲

␮

where qC␮+1i is the position vector of the vertex i of the new curve C␮+1, and the matrix ⵱qCi⵱qT Eapprox共qCi , gCi , ⌬qCi兲兩C=C␮ Ci is the Weierstrass E-Hessian matrix evaluated in the point vertex i of the curve C␮. Using the above projectors and their properties, Eq. 共E3兲 can be rearranged as B共i ␮兲共qC␮+1i − qC␮i兲

关 ⫻兵兩⵱ + 关兩⵱

= − 兩⵱qCi⵱qT Eapprox共qCi,gCi,⌬qCi兲兩 Ci

C=C␮

兴

−1

兩

qCiEapprox共qCi,gCi,⌬qCi兲 C=C ␮

兩


其

兴

⫻A共i ␮兲共qC␮+1i − qC␮i兲 .

共E4兲

From Eq. 共E4兲 we see that if we know the variation of the position vector within the current Krylov subspace, A共i ␮兲共qC␮+1i − qC␮i兲, then we can compute the variation of the position vector in the corresponding complementary subspace B共i ␮兲共qC␮+1i − qC␮i兲. If we define the vector c共i ␮兲 of dimension ␮ − 1 as R共i ␮兲c共i ␮兲 = A共i ␮兲共qC␮+1i − qC␮i兲

共E5兲

and we multiply Eq. 共E4兲 from the left by A共i ␮兲, we obtain

关

0 = A共i ␮兲兩⵱qCi⵱qT Eapprox共qCi,gCi,⌬qCi兲兩

兵 + 关兩⵱

Ci

⫻ 兩⵱qCiEapprox共qCi,gCi,⌬qCi兲兩C=C

C=C␮

兴

−1

␮

兩


兴R

共␮兲共␮兲 i ci

其. 共E6兲

Now, the action of the Weierstrass E-Hessian matrix ⵱qCi⵱qT Eapprox共qCi , gCi , ⌬qCi兲兩C=C␮ on the R共i ␮兲 matrix is apCi proximated as follows:

234105-16

关兩⵱

J. Chem. Phys. 122, 234105 共2005兲


兩


⬇

兴R

关共兩⵱q Eapprox共qCi,gCi,⌬qCi兲兩C=C Ci

共␮兲 i

1

兲

− 兩⵱qCiEapprox共qCi,gCi,⌬qCi兲兩C=C ,…, ␮

共

⫻ 兩⵱qCiEapprox共qCi,gCi,⌬qCi兲兩C=C − 兩⵱qCiEapprox共qCi,gCi,⌬qCi兲兩C=C

␮

␮−1

兲兴 = G共i ␮兲 .

共E7兲

The G共i ␮兲 matrix is a rectangular matrix of dimension N ⫻ 共␮ − 1兲 defined by the set of vector differences between the current and previous Weierstrass E-gradient vectors. Multiplying Eq. 共E6兲 from the left by 共R共i ␮兲兲T and substituting Eq. 共E7兲, we obtain 0

共␮兲

=共

兲关兩

R共i ␮兲 T

兵

兴

兩

⫻ 兩⵱qCiEapprox共qCi,gCi,⌬qCi兲兩C=C + ␮

其, 共E8兲

where 0共␮兲 is a zeroed vector of dimension ␮ − 1. To solve Eq. 共E8兲, first we define the matrix

关

D共i ␮兲 = 共R共i ␮兲兲T 兩⵱qCi⵱qT Eapprox共qCi,gCi,⌬qCi兲兩 Ci

C=C␮

兴

−1

G共i ␮兲共E9兲

and the vector

关

e共i ␮兲 = 共R共i ␮兲兲T 兩⵱qCi⵱qT Eapprox共qCi,gCi,⌬qCi兲兩 Ci

C=C␮

兴

−1

⫻ 兩⵱qCiEapprox共qCi,gCi,⌬qCi兲兩C=C .

共E10兲

␮

Notice that the dimensions of the D共i ␮兲 matrix and the e共i ␮兲 vector are 共␮ − 1兲 ⫻ 共␮ − 1兲 and 共␮ − 1兲 respectively. With the above definitions, the solution of Eq. 共E8兲 is c共i ␮兲 = − 关共D共i ␮兲兲T共D共i ␮兲兲兴−1共D共i ␮兲兲Te共i ␮兲 .

共E11兲

c共i ␮兲

vector in Eq. 共E5兲 Substituting the resulting value of the we obtain A共i ␮兲共qC␮+1i − qC␮i兲 vector. Using the set of equations 共E4兲, 共E5兲, and 共E7兲, we obtain an expression to compute the term B共i ␮兲共qC␮+1i − qC␮i兲, B共i ␮兲共qC␮+1i − qC␮i兲

关 ⫻兵兩⵱

= − 兩⵱qCi⵱qT Eapprox共qCi,gCi,⌬qCi兲兩 Ci

兩

C=C␮

qCiEapprox共qCi,gCi,⌬qCi兲 C=C ␮

兴

−1

其

+ G共i ␮兲c共i ␮兲 . 共E12兲

A共i ␮兲共qC␮+1i − qC␮i兲

and the vector Finally, taking the vector 共␮兲 Bi 共qC␮+1i − qC␮i兲, evaluated from Eq. 共E12兲, we obtain the position vector of the vertex i of the new curve C␮+1, qC␮i + A共i ␮兲共qC␮+1i − qC␮i兲 + B共i ␮兲共qC␮+1i − qC␮i兲 = qC␮i + 共qC␮+1i − qC␮i兲 = qC␮+1i .

W. Quapp, J. Theor. Comput. Chem. 2, 385 共2003兲. J. González, X. Giménez, and J. M. Bofill, J. Phys. Chem. A 105, 5022 共2001兲; Theor. Chem. Acc. 112, 75 共2004兲. 3 K. Fukui, J. Phys. Chem. 74, 4161 共1970兲. 4 W. R. Hamilton, Philos. Trans. R. Soc. London 95, 共1835兲. 5 C. G. J. Jacobi, J. Reine Angew. Math. 17, 97 共1837兲. 6 R. Courant and D. Hilbert, Methods of Mathematical Physics 共Wiley, New York, 1953兲. 7 A. Tachibana and K. Fukui, Theor. Chim. Acta 49, 321 共1978兲. 8 A. Tachibana and K. Fukui, Theor. Chim. Acta 51, 189 共1979兲. 9 A. Tachibana and K. Fukui, Theor. Chim. Acta 51, 275 共1979兲. 10 A. Tachibana and K. Fukui, Theor. Chim. Acta 57, 81 共1980兲. 11 K. Fukui, Int. J. Quantum Chem., Quantum Chem. Symp. 15, 633 共1981兲. 12 L. L. Stachó and M. I. Bán, Theor. Chim. Acta 83, 433 共1992兲. 13 P. G. Mezey, Theor. Chim. Acta 62, 133 共1982兲. 14 A. Ulitsky and R. Elber, J. Chem. Phys. 92, 1510 共1990兲. 15 C. Choi and R. Elber, J. Chem. Phys. 94, 751 共1991兲. 16 R. Olenden and R. Elber, J. Mol. Struct.: THEOCHEM 398-399, 63 共1997兲. 17 H. Jónsson, G. Mills, and K. W. Jacobsen, in Classical and Quantum Dynamics in Condensed Phase Simulations, edited by B. J. Berne, G. Cicotti, and D. F. Coker, 共World Scientific, Singapore, 1998兲, p. 385. 18 G. Henkelman and H. Jonsson, J. Chem. Phys. 113, 9978 共2000兲. 19 W. E, W. Ren, and E. Vanden-Eijnden, Phys. Rev. B 66, 052301 共2002兲. 20 W. Ren, Commun. Math. Sci. 1, 377 共2003兲. 21 S. A. Trygubenko and D. J. Wales, J. Chem. Phys. 120, 2082 共2004兲; 120, 7820 共2004兲. 22 W. Quapp, J. Comput. Chem. 25, 1277 共2004兲. 23 Y. Kanai, A. Tilocca, A. Selloni, and R. Car, J. Chem. Phys. 121, 3359 共2004兲. 24 H. B. Schlegel, J. Comput. Chem. 24, 1514 共2003兲, and references therein. 25 K. Ruedenberg and J.-Q. Sun, J. Chem. Phys. 100, 5836 共1994兲. 26 H. Rund, The Hamilton-Jacobi Theory in the Calculus of Variations 共Van Nostrand, London, 1966兲. 27 W. Yourgrau and S. Mandelstam, Variational Principles in Dynamics and Quantum Theory, 共Dover, New York, 1979兲. 28 B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and M. Karplus, J. Comput. Chem. 4, 187 共1983兲. 29 J.-W. Chu, B.L. Trout, and B. R. Brooks, J. Chem. Phys. 119, 12708 共2003兲. 30 R. Fletcher, Practical Methods of Optimization 共Wiley, New York, 1987兲. 31 M. Page and J. W. McIver, Jr., J. Chem. Phys. 88, 922 共1988兲. 32 C. T. Kelley, Iterative Methods for Linear and Nonlinear Equations, Frontiers in Applied Mathematics 16 共SIAM, Philadelphia, 1995兲. 33 D. R. Fokkema, G. L. G. Sleipjen, and H. A. Van der Vorst, SIAM J. Sci. Comput. 共USA兲 19, 657 共1998兲. 34 R. J. Harrison, J. Comput. Chem. 25, 328 共2004兲. 1 2

−1 ⵱qCi⵱qT Eapprox共qCi,gCi,⌬qCi兲 Ci C=C␮

G共i ␮兲c共i ␮兲

E-Hessian matrix ⵱qCi⵱qT Eapprox共qCi , gCi , ⌬qCi兲兩C=C␮, which Ci appears in Eq. 共E11兲, through the D共i ␮兲 matrix and the e共i ␮兲 vector, and in Eq. 共E12兲, is normally taken as a unit matrix or as a diagonal matrix. As a summary the algorithm at the 共␮ + 1兲th iteration, and for each point vertex i, only needs the position vector q C␮i, the Weierstrass E-gradient vector ⵱qCiEapprox共qCi , gCi , ⌬qCi兲兩C=C␮, and the matrices R共i ␮兲 and G共i ␮兲, defined in the expression 共E2兲 and 共E7兲, respectively. With these vectors and matrices solve the set of equations 共E11兲, 共E5兲, and 共E12兲, and finally using Eq. 共E13兲 compute the new position vector qC␮+1i.

共E13兲

In standard applications of the above method, the Weierstrass