Dynamical dimer method for the determination of transition states with ab initio molecular dynamics
aa r X i v : . [ phy s i c s . c o m p - ph ] N ov Dynamical dimer method for the determination of transition states with ab initiomolecular dynamics
Alexander Poddey ∗ and Peter E. Bl¨ochl † Clausthal University of Technology, Institute of Theoretical PhysicsLeibnizstrasse 10, D-38678 Clausthal-Zellerfeld, Germany (Dated: October 25, 2018)A dynamical formulation of the dimer method for the determination of transition states is pre-sented. The method is suited for ab-initio molecular dynamics using the fictitious Lagrangian for-mulation. The method has been applied to the con-rotatory ring opening of chloro-cyclo-butadiene,an example, where the application of the drag method is problematic.
I. INTRODUCTION
The concept of the transition state has a fundamen-tal role for the prediction of rate constants of materi-als processes. The transition state is the lowest pointon the energy barrier separating two metastable statesrepresenting the initial and final state of a process. Theenergy difference between the initial state and the transi-tion state is the activation energy. The curvatures of thepotential energy surface at the initial state and the tran-sition state provides an estimate of entropic effects[1, 2].A molecular dynamics simulation starting from the tran-sition state allows one to construct the dynamical tra-jectory of a reaction in the low-temperature limit. Theextension of this latter approach to finite temperaturesis straightforward [3].The determination of transition states is however sub-stantially more difficult than the determination of min-ima of the total energy surface. A large number of meth-ods have been invented and refined over the last decades(see e.g. the Refs. 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15and the review articles [16, 17]).Many of these methods require the determination ofthe second derivatives of the total energy surface. Whileforces, that is first derivatives of the potential energy sur-face, are readily accessible from density functional calcu-lations, second derivatives are not, at least not withoutsubstantial effort.An overview of methods which allow to determine tran-sition states solely from force information can be foundin [16]. The most simple approach is the drag method,where a one-dimensional constraint forces the systemacross a reaction barrier. Widely used are also “chainof states” algorithms such as the nudged elastic band(NEB) [10, 11, 12] or the string method [14, 15].Intermediate between these extremes is the dimermethod. [18]. It couples two instances of the system,which we will call monomers. The monomers have aspecified distance in configuration space. If the systemfollows the forces, after inverting of the force parallel to ∗ Electronic address: [email protected] † Electronic address: [email protected] the dimer, it evolves to a saddle point of first order, thatis the transition state.The original dimer method [18] applied a two stepalgorithm. In a first step, the potential energy of thedimer-construct is minimized, excluding parallel motion,following a conjugate gradient approach. Then a largestep is taken in the direction parallel to the dimer axis.Recently, an improvement has been suggested which min-imizes the number of gradient evaluations [19].In the context of fictitious-Lagrangian approach tofirst-principles molecular dynamics (FPMD) [20], such atwo-step procedure is not feasible, because the electronicwave functions must be able to follow the atomic mo-tion adiabatically. The fictitious Lagrangian approachrequires a larger number of time steps. However thecomputational cost per time step are minimal. Whilea first-principles string approach has been published re-cently [21], the dimer method has not been examined inthe framework of ab initio molecular dynamics.In the present study, we introduce a dynamical formu-lation of the dimer method. Starting from an extendedLagrangian for the dimer in the configuration space withdoubled dimensionality, we derive equations of motion,that, in the presence of dissipation, converge at transi-tion states of first order. An analysis of the trajectoriesof the dimer close to the stable points provides informa-tion on the stability and gives guidance for improving theperformance of the method. The method has been im-plemented in the CP-PAW code and has been applied tothe con-rotatory ring-opening of chloro-cyclo-butene, anexample where the conventional drag method fails unlessspecial precautions are taken.The paper is organized as follows: In section II we de-fine the extended Lagrangian for the dimer motion andderive the equations of motion. In section III we ana-lyze stability and the dynamics in the proximity of thestable points. Section IV is devoted to the discretizationthe equations of motion and the aspects of the actualimplementation. In the final section V, the method isapplied to the con-rotatory ring opening of cyclo-buteneas practical example and compared to the drag method.
II. LAGRANGIAN AND EQUATIONS OFMOTION
In this section we define the Lagrangian underlyingthe present work and derive the equations of motion. Webegin by defining our notation and we introduce a con-venient set of coordinates:The dimer consists of two points in configurationalspace separated by a fixed distance. We use mass-weighted coordinates x i defined as x i := m R i , (1)where the vectors R i with i = 1 , N atoms each vector has 3 N dimensions. In addition, each monomer has its own set ofelectronic wave functions. The mass matrix m is diagonaland the masses of the respective atoms are located on themain diagonal.It is convenient to introduce a variable transformationinto center-of-gravity coordinates Y and relative coor-dinates Y . Y = x + x Y = x − x (3)Loosely speaking Y describes the mean structure of themolecule, while Y describes the orientation of the dimer.A dimer has three basic types of motion: motion par-allel to the dimer axis ( v k ), motion perpendicular to theaxis ( v ⊥ ) and rotation about the center of gravity ( v ◦ ).In the following, we divide the velocities into the con-tributions from these basic types of motion. The corre-sponding velocities are v ◦ = ˙Y = ˙x − ˙x (4) v k = Y ⊗ Y d ˙Y = ( x − x ) ( x − x ) ( ˙x + ˙x )2 d (5) v ⊥ = (cid:18) − Y ⊗ Y d (cid:19) ˙Y = ˙x + ˙x − v k (6)The dot denotes the time derivative, and d is the so called’dimer distance’ d = p Y = q ( x − x ) . The sym-bol ⊗ refers to the outer product defined by ( a ⊗ b ) c = a ( b · c ). Note, that Y ⊗ Y d is an operator that projectsa vector onto the dimer axis.Now we define our theory by setting up a Lagrangianfunction in the relative and center-of-gravity coordinates L ( Y , Y , ˙Y , ˙Y ) = M ◦ v ◦ + M ⊥ v ⊥ + M k v k − V (cid:18) m − (cid:20) Y + 12 Y (cid:21)(cid:19) − V (cid:18) m − (cid:20) Y − Y (cid:21)(cid:19) − ¯ λ (cid:2) Y − d (cid:3) (7) The potential energy is described by the potential en-ergy surface V ( R ). The velocities v ◦ , v ⊥ and v k areto be treated as functions of Y and Y and their timederivatives ˙Y and ˙Y as defined in Eqs. 4,5,6.We introduced three scale factors M ◦ , M k and M ⊥ ,that allow us to scale the masses for the three types ofmotions independently, that will later be used to accel-erate convergence. More importantly, by choosing a neg-ative value of M k , the motion along the dimer directionis inverted, so that the dimer climbs up the barrier.If the scale factors are equal to one, the conventionalkinetic energy is recovered, that is14 v ◦ + v k + v ⊥ = 12 ˙R m ˙R + 12 ˙R m ˙R (8)Using the method of Lagrange multipliers, a termhas been introduced to describe the dimer-distance con-straint Y = d (9)which ensures that the dimer distance in mass-weightedcoordinates is equal to d . The corresponding Lagrangemultiplier is denoted by ¯ λ .In addition to the Lagrangian we also define aRayleigh’s dissipation function DD = γ ◦ M ◦ v ◦ + γ k M k v k + γ ⊥ M ⊥ v ⊥ . (10)which allows us to introduce dissipation in a consistentmanner.The Euler-Lagrange equation for the Lagrangian ofEq. 7 and the Rayleigh’s Dissipation function from Eq. 10are obtained from ddt ∂ L ∂ ˙ Y i,n − ∂ L ∂Y i,n + ∂ D ∂ ˙ Y i,n = 0 , (11)where the index i ∈ { , } labels the two monomers, and n ∈ { , N } labels the coordinate in configuration space.The resulting equations of motion have the form (cid:18) − Y ⊗ Y d (cid:19) M ⊥ (cid:16) ¨ Y + γ ⊥ ˙Y (cid:17) + (cid:18) Y ⊗ Y d (cid:19) M k (cid:16) ¨ Y + γ k ˙Y (cid:17) − m − ( F + F )+ (cid:0) M k − M ⊥ (cid:1) (cid:20) ddt (cid:18) Y ⊗ Y d (cid:19)(cid:21) ˙Y = 0 (12)and M ◦ ¨ Y − m − ( F − F ) + 4 λ Y +4 (cid:0) M ⊥ − M k (cid:1) Y ˙Y d ˙Y + γ ◦ M ◦ ˙Y = 0 (13)Here we have used the forces acting on the monomersdefined as F i := − ∇| R i V (14)In equation (13) we absorbed all terms which lead to aforce parallel to the dimer axis into the constraint forceby redefining the Lagrange multiplier. Thus the variable λ used in Eq. 13 differs from the Lagrange multiplier ¯ λ in Eq. 7.Furthermore, because the dimer system moves on thehyperplane with constant dimer distance, we simplifiedthe final equations by using Y = d (15)and ∂∂t Y = 0 . (16)Equation (12) describes the motion of the center of grav-ity of the dimer, which is related to the structure of themolecule. Equation (13) describes the orientational mo-tion of the dimer. III. LOCAL STABILITY ANALYSIS
For M ◦ = M k = M ⊥ = 1, the Lagrange function (7)describes a physical system of two masspoints moving ina potential V under the influence of the constraint forcethat keeps the dimer distance invariant. With positivefriction factors γ ◦ , γ k and γ ⊥ the dimer will come to restwith the center-of-gravity coordinate next to a local min-imum. The dimer axis will be aligned nearly parallel tothe lowest vibrational eigenmode. The above statementsare exactly fulfilled in the limit of vanishing dimer dis-tance.If we choose a negative value of M k , the motion will be-come unstable near local minima. Instead, the dimer willbe attracted by transition states of first order. A transi-tion state of first order is characterized by the presenceof exactly one eigenmode with an imaginary frequency.These properties of the dynamics have been derivedfrom the following local stability analysis. First we de-termine the stationary points for the dimer dynamics.Then we investigate the dynamics in the neighborhoodof those stationary points.For this analysis we can replace the potential energysurface by its truncated Taylor expansion at a given point R (0) = m − x (0) , namely V ( R ) = V (0) − F (0) R + 12 ( R − R (0) ) m Dm ( R − R (0) ) (17)with the dynamical matrix D defined as D i,j = 1 √ m i,i ∂ V∂R i ∂R j (cid:12)(cid:12)(cid:12)(cid:12) R (0) √ m j,j (18) From Eq. 17 and Eq. 1, we obtain the forces F ( x ) = F (0) − m D ( x − x (0) ) (19) A. Stationary points
If we insert the condition for a stationary point ˙Y i = ¨Y i = 0 in the equation of motion Eq. 12, we find thatthe forces for the two monomers at the stationary pointare antiparallel and of the same magnitude. Using theTaylor expansion Eq. 19 in the equation of motion Eq. 12,we find that the center-of-gravity of the dimer Y liesat a stationary point of the potential, if the dimer isstationary, that is F ( Y ) = 0 (20)Similarly we obtain from Eq. 13 and Eq. 19 DY = − λ Y (21)Thus the dimer axis points along an eigenvector of thedynamical matrix.In conclusion we find that the dimer dynamics is sta-tionary, when its center of gravity lies at an extremum ora saddle point of the potential energy surface and when,in addition, the dimer axis points along one of the vibra-tional eigenmodes. B. Linearized equation of motion
Without loss of generality, we choose in the followinga coordinate system for which the stationary point of thepotential energy surface lies at the origin. Secondly, wedenote the eigenvectors of the dynamical matrix as e i ,with e i e j = δ i,j , and the corresponding eigenvalues with ω i . The eigenvector, which is parallel to the dimer axis isdenoted by ¯e and the corresponding eigenvalue is denotedby ¯ ω . Thus, with Eq. 21, we can identify the Lagrangemultiplier as λ = − ¯ ω .Linearization of the equations of motion Eqs. 12 and13 about Y = and Y = ¯e d with F (0) = 0 yields thefollowing equations for the deviation δ Y ( t ) and δ Y ( t )from the stationary point(1 − ¯e ⊗ ¯e ) M ⊥ (cid:16) δ ¨ Y + γ ⊥ δ ˙ Y (cid:17) + ( ¯e ⊗ ¯e ) M k (cid:16) δ ¨ Y + γ k δ ˙ Y (cid:17) + D δ Y = 0 (22) M ◦ δ ¨ Y + D δ Y + ¯ ω δ Y + M ◦ γ ◦ δ ˙ Y = 0 (23)We project the second equation ,Eq. 23, onto the eigen-vectors e i and obtain: M ◦ ( e i δ ¨ Y ) + ( ω i − ¯ ω )( e i δ Y ) + M ◦ γ ◦ ( e i δ ˙ Y ) = 0(24)Note, that the motion along ¯e is simple due to the dis-tance constraint. The frequency of the variable e i Y istherefore ω ◦ = ± s ω i − ¯ ω M ◦ (25)We see that, for a positive mass M ◦ , the dynamics isstable only, if ¯ ω is the lowest eigenvalue of the dynamicalmatrix. Thus the dimer axis will always orient along theeigenvector with the lowest eigenvalue.Now, we project Eq. 22 onto the eigenvector ¯e of thedynamical matrix M k ( ¯e δ ¨ Y ) + ¯ ω ( ¯e δ Y ) + M k γ k ( ¯e δ ˙ Y ) = 0 (26)which yields the translational motion of the dimer alongthe dimer axis. This equation will be responsible forthe ascent to a saddle point. The eigenfrequency of thevariable ¯eY is ω k = ± ¯ ω p M k (27)For a negative mass M k , which is the choice for a tran-sition state search, Eq. 26 is stable only if the dimer isoriented along an unstable vibrational mode of the po-tential energy surface, that is with the dimer orientedalong an eigenvector ¯e with an imaginary frequency.Now, we project Eq. 22 onto the eigenvectors e i with e i ⊥ ¯e and we obtain M ⊥ ( e i δ ¨ Y ) + ω i ( e i δ Y ) + M ⊥ γ ⊥ ( e i δ ˙ Y ) = 0 (28)which represents the translational motion perpendicularto the dimer axis. The translational motion is related toan optimization of the atomic structure. The frequencyof the variable e i Y with e i ⊥ ¯e is ω ⊥ = ω i √ M ⊥ (29)From Eqs. 25, 27 and 29 we find that the dimer withpositive M ◦ and M ⊥ and negative M k will only come torest at transition states of first order, that is, if ¯ ω < C. Masses and Frictions
We will later solve the equations of motion 12, 13 withthe Verlet algorithm. The Verlet algorithm for a har-monic oscillator with a frequency ω becomes unstable, ifthe time step ∆ used for the discretization of the equa-tion of motion is smaller than 2 /ω . The frequency of thediscretized motion is accurate to within 1 % if ω ∆ < π .Thus, it is important to understand the vibrational spec-trum, shown schematically in Fig. 1, of the motion in thepotential. ω ω ω min2 ω max2 (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) FIG. 1: Sketch of the eigenvalue spectrum of the dynamicalmatrix at the saddle point of first order with the definition of ω min and ω max . Thus, the stability limit requires us to choose themasses as M k < ¯ ω ∆ M ⊥ > ω max ∆ M ◦ > ( ω max − ω min )∆ a priori estimate of these lim-its: (1) For practical purposes, we can make the assump-tion that ω min << ω max and thus we set ω min to zero.(2) The highest vibrational mode is that of H , with afrequency of 4400 cm − . Frequencies for bond-stretchvibrations are of the order 1000 cm − . These numbersare reasonable upper estimates for ω max . (3) The abso-lute value of the imaginary frequency at the saddle pointis typically in the range of the high-lying real-frequencymodes, that is comparable to ω max . Hence ¯ ω ≈ − ω max .With these assumptions we obtain the following, recom-mended values for the masses M k = − κ ∆ (33) M ⊥ = κ ∆ (34) M ◦ = 2 κ ∆ (35)where κ = (cid:0) · − a . u . (cid:1) is a recommended constantbased on the H vibration.Let us now turn our attention to the friction values.The goal is to reach an optimum convergence at the sta-tionary point of the dimer. For a damped harmonic os-cillator m ¨ x = − mω x − mγx (36)the fastest decay rate as function of the friction is ob-tained for critical damping, that is in between the oscil-latory and the over-damped regime. Critical damping isobtained for a friction of γ = 2 ω . Thus we choose γ k = 2 s ¯ ω M k (37)2 s ω min M ⊥ < γ ⊥ < s ω max M ⊥ (38)2 s ω min − ¯ ω M ◦ < γ ◦ < s ω max − ¯ ω M ◦ (39)Note, that a too high friction freezes out those modes inthe over-damped regime. Thus it is usually better to usea friction near the lower bound of the sensible regime. IV. NUMERICAL INTEGRATION OF THEEQUATIONS OF MOTIONA. Discretization of the equations of motion
The equations of motion Eqs. 12,13 are nonlinear inthe velocities. This leads to a nonlinear equation for thediscretized equation of motion. We tackle the problemby iterating on the nonlinear terms in the velocities. We set up the discretized equations of motion, whiletreating the terms nonlinear in the velocities as an ab-stract force. G := − (cid:0) M k − M ⊥ (cid:1) (cid:20) ddt (cid:18) Y ⊗ Y d (cid:19)(cid:21) ˙Y (40) G := 4 (cid:0) M k − M ⊥ (cid:1) Y ˙ Y d ˙Y (41)The equations of motion are discretized according tothe Verlet algorithm:˙ y → y (+) − y ( − )2∆ (42)¨ y → y (+) − y (0) + y ( − )∆ (43)where ∆ is the discretization time step, and where weused the short-hand notation y (+) = y ( t + ∆) (44) dy (0) = y ( t ) (45) y ( − ) = y ( t − ∆) (46)We obtain from Eq. 12 (cid:18) − Y ⊗ Y d (cid:19) M ⊥ ∆ n (1 + a ⊥ ) Y (+) − Y (0) + (1 − a ⊥ ) Y ( − ) o + (cid:18) Y ⊗ Y d (cid:19) M k ∆ n(cid:0) a k (cid:1) Y (+) − Y (0) + (cid:0) − a k (cid:1) Y ( − ) o − m − ( F + F ) − G = 0 (47)where a ⊥ = γ ⊥ ∆2 and a k = γ k ∆2 . In the projectors wehave dropped the argument of Y (0).Eq. 47 can be resolved for Y (+) by multiplicationwith (cid:18) − Y ⊗ Y d (cid:19) ∆ M ⊥ (1 + a ⊥ )+ (cid:18) Y ⊗ Y d (cid:19) ∆ M k (1 + a k ) (48) The result is Y (+) = (cid:18) − Y ⊗ Y d (cid:19) n
21 + a ⊥ Y (0) − − a ⊥ a ⊥ Y ( − ) + (cid:18) m − ( F + F ) + G (cid:19) ∆ M ⊥ (1 + a ⊥ ) o + (cid:18) Y ⊗ Y d (cid:19) n
21 + a k Y (0) − − a k a k Y ( − ) + (cid:18) m − ( F + F ) + G (cid:19) ∆ M k (1 + a k ) o (49)Similarly we discretize the equation of motion for the dimer distance, Eq. 13. Following Ryckaert et al [22], wefirst propagate without the force of constraint to obtain¯ Y = 21 + a ◦ Y (0) − − a ◦ a ◦ Y ( − )+ (cid:16) m − ( F − F ) + G (cid:17) ∆ M ◦ (1 + a ◦ ) (50)with a ◦ = γ ◦ ∆2 . The new vector Y (+) is related to ¯ Y by the constraint force and the Lagrange parameter λ according to Y (+) = ¯ Y − λ Y (0) ∆ M ◦ (1 + a ◦ ) (51)The Lagrange parameter is then adjusted so that theconstraint Y (+) = d , namely a given dimer length, issatisfied.In the first iteration we estimate this forces G and G from the previous iterations or we set them to zero.Then we propagate the positions, which provides us witha better estimate for the nonlinear terms. This loop isthen iterated to convergence. B. Restricting the orientation of the dimer
It will be important to limit the degrees of freedom,in which the two monomers may differ. This implies re-stricting the orientation of the dimer to a space withlower dimensionality. The reason is to avoid that thedimer converges at irrelevant saddle-points, which are notof interest. This is a common problem for complex sys-tems, that exhibit many minima and saddle points. Thisrestriction is accomplished easily, by setting the corre-sponding components of Y in Eq. 50 to zero. V. APPLICATION: CONROTATORY RINGOPENING OF 1-CHLORO-2-CYCLOBUTENE
In order to demonstrate the performance of the dimermethod described above, we have chosen a prototypicalsystem for tests of transition state searches, namely thering opening of cyclobutene. The isomerization of cy-clobutene to cis-butadiene is the prototypical example ofconcerted stereospecific reactions. The underlying pro-cesses have been studied and discussed extensively byseveral groups [23, 24, 25, 26] (and references therein).The potential energy surfaces of the con- as well as ofthe disrotatory isomerization mechanisms exhibit prin-cipal structures, which render the determination of thecorresponding transition states complicated [24].In order to avoid special effects due to the high sym-metry of cyclobutene, that would be untypical for othermolecules, we explored 1-chloro-2-cyclobutene. Thismolecule and the two relevant reaction products areshown in Figure 2.Computational details of our calculations are given inappendix A. The coordinates of the structures are sup-plied with the supplementary material.
Cl C A CB C C C C H H H H H H C C H H C H H C H C C H C H H H ClCl
FIG. 2: Initial state, 1-chloro-2-cyclobutene (A), and finalstate, 1-chloro-buta-1,3-diene (B), for the conrotatory ringopening of chlorocyclobutene. Also shown is the productof the disrotatory ring opening, trans-1-chloro-buta-1,3-diene(C).
A. Drag Calculations
In this section we investigate the reaction using thedrag method. Because of its simplicity, the drag methodis widely used for transition state search. However, thedrag method fails for a number of systems. In such casesthe energy exhibits a hysteresis effect, that is, differentenergy profiles are obtained when dragging the reactioncoordinate in the opposite directions. The ring openingof cyclobutene is one example where the drag methodexhibits a hysteresis effect.In the drag method one specifies a one-dimensionalreaction coordinate. Then one maps the energy profilealong the reaction coordinate and determines the highestpoint as the transition state. Each point along the en-ergy profile corresponds to the local energy minimum ona hyperplane perpendicular to the reaction coordinate.In our case we selected the reaction coordinate by thedifference between initial and final state. A hyperplanewith a specified value c of the reaction coordinate is givenby ( R B − R A ) ( R − R A )( R B − R A ) = c (52)Hence the value of the reaction coordinate is zero for theinitial state (A) and one for the final state (B).In addition to the reaction coordinate we imposed sixadditional constraints to avoid translations and rotationsof the molecule. We have chosen12 ( R ( C ) + R ( C )) = 0 (53) z ( C ) = z ( C ) = 12 ( x ( C ) + x ( C )) = 0 (54)Where the coordinates correspond to the atoms given inparenthesis. The notation follows Figure 2.When the energy profile is determined by varying thereaction coordinate in small steps, once from zero to oneand then from one to zero, we obtain the hysteresis shownin Fig. 3. The highest point of the energy profile is notthe transition state. Instead, the highest point of eachbranch is an upper bound for the activation energy, whilethe crossing of the two branches is a lower bound. E ne r g y ( e V ) FIG. 3: Potential energy relative to the energy of final state(B) over reaction coordinate from a forward and backwardcalculation using the drag method. The true transition state(determined using the dynamical dimer method) is denotedby a circle.
In order to show the underlying reason for the hys-teresis, we show in Fig. 4 the total energy surface in atwo-dimensional hypersurface. The two axes are the re-action coordinate and the coordinate defined by the twoisoenergetic structures at the crossing of the two branchesof the energy profile in Fig. 3. For each point in the two-dimensional plot the energy is at a local minimum of the(3 N − A † and B † , that markthe crossing of the two branches in the energy profileshown in Fig. 3, lie far from the transition state. Forfurther information see Table I. AA c c C BB TS FIG. 4: Sketch of the potential energy surface of the con-rotatory isomerization of 1-chloro-2-cyclo-butene to 1-chloro-buta-1,3-diene. The dashed lines trace the paths obtainedfrom drag calculations for the transitions A → B and B → A,respectively. The crosses denote the structures A † and B † .TS corresponds to the transition state. B. Dynamical dimer calculations
In this section we describe technical issues for dy-namical dimer calculations and demonstrate their per-formance.We use the fictitious Lagrangian formulation of ab-initio molecular dynamics. This implies that wave func-tion coefficients and nuclei obey Newton’s equation ofmotion, to which we added a friction, which allows toquench the system.Special attention should be given to the wave functiondynamics. The wave-function cloud tied to the atoms re-sults in an increased effective mass of the nuclei. Further-more, a friction applied to the wave function dynamicsacts like an effective friction on the nuclear motion. Fora normal ground-state search, this is not problematic.In our case however, the mass for the motion parallel tothe dimer direction is inverted. Thus the dimer acceler-ates opposite to the direction of the force. Thus a frictionacting on the dimer via the electrons leads to an velocity-dependent acceleration of the dimer, that may cause aninstability of the dimer motion.The detrimental effect of the wave function motion canbe avoided by choosing a sufficiently small mass for thewave function dynamics or by artificially increasing theatomic masses. In our study, we used a mass of 50 u forall atoms. The wave functions have been kept close tothe ground state by applying a friction that dissipates2 % of the kinetic energy in each time step.The dimension-less masses for the dimer motion havebeen chosen to M k = − .
00 (55) M ⊥ = +1 .
00 (56) M ◦ = +0 . . (57)The small value for M ◦ has been chosen to speed up thereorientation of the dimer.In order to avoid instabilities, we found it useful tointroduce an upper limit to the kinetic energy individu-ally for the rotational, perpendicular and parallel motion.These limits were enforced by adjusting the frictions ac-cordingly. The limit for the kinetic energy E kin,max istranslated into an upper “temperature” T max accordingto 12 gk B T max = E kin,max (58)where g is the number of degrees of freedom in the cor-responding motion type, k B is Boltzmann’s constant.Within the limits of the enforced maximum kinetic en-ergies, we adapted the friction dynamically in each stepto come close to the limit of critical damping. The latterresults in the best possible convergence behavior.Furthermore, the friction is increased to a higher valuespecified by the user, if the dimer moves away from thestable point. The decision to increase the friction hasbeen made on the basis of the actual forces and velocities.In our simulation we have chosen this friction so that theenergy dissipated per time step corresponds to 20 % ofthe kinetic energy of the corresponding type of motion.When starting the dimer simulation, it is importantto optimize the dimer orientation, i.e. Y , first, beforethe mean dimer configuration Y changes appreciably.This is because the mean dimer configuration only ap-proaches the transition state, if the dimer orientation issufficiently aligned along the unstable vibrational modeof the transition state.The dimer orientation can be optimized in differentways. (1) Starting from initial dimer coordinates, themean dimer configuration, Y , is constrained to the ini-tial value, while the orientation, Y , is fully optimized.Only after this optimization also the mean configurationis allowed to move. (2) Starting from one monomer con-figuration, the second monomer is constructed from thefirst by letting it follow the forces until the desired dimerlength is obtained. We will refer to this technique asthe “growing-dimer technique” (3) Starting from a dimerwith the monomers being identical to the two metastableminima, the dimer length is successively shortened, whilethe dimer position is optimized. Results for all three op-timization strategies can be found in Table I.We found that the third strategy suffers from the factthat additional effort is required to contract the dimerlength to the desired value. It also suffers from an insta-bility caused by strong anharmonic effects. The growing-dimer technique has the advantage thatthe second monomer is constructed dynamically from thefirst, without the need of an independent optimization ofthe wave functions. In the growing dimer technique thedimer orientation is optimized automatically as soon asthe targeted dimer length is reached. In the first strategythe orientation is obtained during the first iterations witha constrained mean dimer configuration.We will discuss here the first strategy. The initialstructure of the dimer is given by a mean dimer config-uration midway between the two metastable states (A)and (B). The dimer length is chosen such that the dis-tance of the two monomers in coordinate space is 0.33 ˚A.During the first 700 time steps, the dimer orientationhas been optimized. Here the mean dimer configurationhas not been constrained, but the maximum temperaturefor the perpendicular motion has been kept at a smallvalue of 10 K. For the parallel and rotational motion,the maximum temperature has been set to 500 K. Afterthe first 700 time steps the maximum temperature forthe perpendicular motion has been increased to 500 K.In Figs. 5,6 we show the convergence of the forces forthe parallel, rotational and perpendicular motion. Goodconvergence is reached after 2000 time steps. The esti-mates for the transition state energy and the deviation ofthe mean dimer configuration from the transition stateis given in Table I. F o r c e ( a . u . ) FIG. 5: Forces over time-step for a calculation following op-timzation strategy (1) described in the text. The dashed anddash-dotted line show the forces parallel to the dimer axis forimage one and two, respectively. The full line correspondsto the resulting force acting on the center of gravity of thedimer. The logarithmic scale of the inset provides detailedinformation for the latter.
VI. SUMMARY
A formulation of the dimer method for searching tran-sition states is presented which can be used in ab-initio F o r c e ( a . u . ) FIG. 6: Rotational (full line) and perpendicular (dashed line)part of the forces acting on the dimer over time-steps for acalculation following optimzation strategy (1) described in thetext.TABLE I: Predicted transition-state energy (relative to finalstate (B)) and deviation from the exact transition state forthe con-rotatory isomerization of 1-chloro-2-cyclo-butene to1-chloro-buta-1,3-diene. The results include three differentdynamical dimer optimization strategies ((1), (2) and (3))described in the text and the forward as well as the backwarddrag calculation shown in Fig. 3.Predicted Energy (eV) ˛˛ R TS − R ˛˛ (˚A)Transition State a A → B B → A a Determined by the use of a well converged dynamical dimer cal-culation with final dimer distance of 0.125 ˚A. molecular dynamics simulations using a fictitious La-grangian. The dimer method is successful in cases wherethe widely used drag method fails. We investigate thedynamics close to the stable points of the dimer analyti-cally. In addition we demonstrate the performance of ourimplementation using a practical example typically usedas test case for transition state search algorithms, namelythe con-rotatory ring opening of (chloro-) cyclo-butene.
APPENDIX A: COMPUTATIONAL DETAILS
We performed density-functional calculations [27,28] based on the projector augmented wave (PAW)method[29, 30]. The gradient-corrected PBE[31] func-tional was used for exchange and correlation. The PAWmethod is a frozen-core all-electron method. Like otherplane-wave based methods, the PAW method leads tothe occurrence of artificial periodic images of the struc-tures. This effect was avoided by explicit subtractionof the electrostatic interaction between them.[32] Wavefunction overlap was avoided by choosing the unit celllarge enough to keep a distance of more than 6 ˚A be-tween atoms belonging to different periodic images. Weused a plane wave cutoff of 30 Ry for the auxiliary wavefunctions of the PAW method. The following sets of pro-jector functions were employed, Cl 2s2p1d, C 2s2p1d, H2s1p, which provides the number of projector functionsper angular momentum magnetic quantum number m ineach main angular momentum channel ℓ .Atomic structures were optimized by damped Car-Parrinello[33] molecular dynamics. We used a time-stepof 10 a.u. (2.5 fs) for all except the dynamical dimer cal-culations. See section V B for further details. The con-vergence was tested by monitoring if the total energychange remains below 10 − Hartree during a simulationof 500 time steps. During the simulation for the conver-gence test, no friction was applied to the atomic motionand the friction on the wave function dynamics was cho-sen sufficiently low to avoid a noticeable effect on theatomic motion.0 [1] H. Eyring, J. Chem. Phys. , 107 (1934).[2] G. Vineyard, J. Phys. Chem. Solids , 121 (1957).[3] J. Keck, Discuss. Faraday Soc. , 1962 (1962).[4] D. Poppinger, Chem. Phys. Lett. , 550 (1975).[5] C. J. Cerjan and W. H. Miller, J. Chem. Phys. , 2800(1981).[6] J. Baker, J. Comp. Chem. , 385 (1986).[7] J. Simons, P. Jorgensen, H. Taylor, and J. Ozment, J.Phys. Chem. , 2745 (1983).[8] A. Banerjee, N. Adams, J. Simons, and R. Shepard, J.Phys. Chem. , 52 (1985).[9] L. J. Munro and D. J. Wales, Phys. Rev. B , 3969(1999).[10] H. Jonsson, G. Mills, and K. W. Jacobsen, Classical andQuantum Dynamics in Condensed Phase Simulations ,Ed. B.J. Berne, G.Ciccotti and D.F.Coker, World Sci-entific p. 385 (1998).[11] G. Mills and H. Jonsson, Phys. Rev. Lett. , 1124(1994).[12] G. Mills, H.Jonsson, and G. Schenter, Surf. Sci. , 305(1995).[13] R. Malek and N. Mousseau, Phys. Rev. E , 7723(2000).[14] W. E, W. Ren, and E. Vanden-Eijnden, Phys. Rev. B ,52301 (2002).[15] B. Peters, A. Heyden, A. T. Bell, and A. Chakraborty,J. Chem. Phys. , 7877 (2004).[16] G. Henkelman, G. Johannesson, and H. Jonsson, Progresson Theoretical Chemistry and Physics , Ed. S. D.Schwartz, Kluwer Academic Publishers p. 269 (2000).[17] H. B. Schlegel, J. Comput. Chem. , 1514 (2003). [18] G. Henkelman and H. Jonsson, J. Chem. Phys. , 7010(1999).[19] R. A. Olsen, G. H. Kroes, G. Henkelman, A. Arnoldsson,and H.Jonsson, J. Chem. Phys. , 9776 (2004).[20] R. Car and M. Parrinello, Phys. Rev. Lett. , 2471(1985).[21] Y. Kanai, A. Tilocca, A. Selloni, and R. Car, J. Chem.Phys. , 3359 (2004).[22] J.-P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, J.Comp. Phys. , 327 (1977).[23] R. B. Woodward and R. Hoffmann, J. Am. Chem. Soc. , 395 (1965).[24] M. J. S. Dewar and S. Kirschner, J. Am. Chem. Soc. , 4292 (1971).[25] J. Breulet and H. F. Schaefer, J. Am. Chem. Soc.(1981).[26] J. M. Olivia, J. Gerratt, P. B. Karadakov, and D. L.Cooper, J. Chem. Phys. , 8917 (1997).[27] P. Hohenberg and W. Kohn, Phys. Rev. , B864(1964).[28] W. Kohn and L. Sham, Phys. Rev. , A1133 (1965).[29] P. E. Bl¨ochl, Phys. Rev. B , 17953 (1994).[30] P. E. Bl¨ochl, C. F¨orst, and J. Schimpl, Bull. Mater. Sci. , 33 (2003).[31] J. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. , 3865 (1996).[32] P. E. Bl¨ochl, J. Chem. Phys. , 7422 (1986).[33] R. Car and M. Parrinello, Phys. Rev. Lett.55