Efficient Transition State Optimization of Periodic Structures through

Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

pubs.acs.org/JCTC

Efficient Transition State Optimization of Periodic Structures through Automated Relaxed Potential Energy Surface Scans Philipp N. Plessow* Institute of Catalysis Research and Technology (IKFT), Hermann-von-Helmholtz-Platz 1, D-76344 Eggenstein-Leopoldshafen, Germany S Supporting Information *

ABSTRACT: This work explores how constrained linear combinations of bond lengths can be used to optimize transition states in periodic structures. Scanning of constrained coordinates is a standard approach for molecular codes with localized basis functions, where a full set of internal coordinates is used for optimization. Common plane wavecodes for periodic boundary conditions almost exlusively rely on Cartesian coordinates. An implementation of constrained linear combinations of bond lengths with Cartesian coordinates is described. Along with an optimization of the value of the constrained coordinate toward the transition states, this allows transition optimization within a single calculation. The approach is suitable for transition states that can be well described in terms of broken and formed bonds. In particular, the implementation is shown to be effective and efficient in the optimization of transition states in zeolite-catalyzed reactions, which have high relevance in industrial processes.

■

others that typically also employ internal coordinates.15,16 Additionally, with the availability of internal coordinates one can define linear combinations of these, for example of bond lengths. This allows the definition of internal coordinates, also called distinguished coordinates, that resemble a reaction coordinate.16−24 One can now perform energy minimizations, while constraining these internal coordinates. Performing constraint optimizations for several values of the constraint allows to “scan” along the approximate reaction coordinate and to thereby obtain the transition state. While this approach has been used in molecular quantum chemistry for a long time, it has never really caught on in surface science or heterogeneous catalysis. The main exception here is the breaking of single bonds on the metal surface, where the constrained variation of this single bond has been used, sometimes referred to as the “fixed bond-length approximation”.25−27 Geometry optimizations with fixed bond lengths can be carried out with the atomistic simulation environment, ASE,12 or VASP. The use of internal coordinates for periodic structures has also been explored.28−30 Recently, Jafari et al. reported an extension to the growing string method (GSM)31 in which internal coordinates are introduced. If a reaction path is now grown from one side, using an initial direction specified in internal coordinates, this approach can be used to drive reactions similar to relaxed potential energy surface scans. However, to the best of my knowledge, linear combinations of bond lengths are typically not used to optimize transitions states in PBC-codes.

INTRODUCTION The calculation of rate constants based on ab initio and density functional theory (DFT) calculations, for example in computational catalysis, largely relies on harmonic transition state theory. This requires the optimization of first-order saddle points that connect reactant and products via a minimum energy path (MEP). Transition state optimization is generally much more difficult than energy minimization.1−3 This challenge along with the demand for efficient algorithms continues to motivate researchers to develop new approaches. Traditionally, optimization algorithms have to some extent always been tied to specific codes or program packages, such as Gaussian,4 Turbomole,5 Qchem,6 and GAMESS-US.7 Furthermore, development of algorithms for the above-mentioned codes that employ atom-centered localized basis functions has followed different directions from those in programs for periodic boundary conditions (PBC) that usually employ plane waves, such as VASP,8,9 Quantum Espresso,10 and GPAW.11,12 This may to some extent be due to little overlap between the communities in which these programs are used. Another factor is the intrinsic difference in the studied systems, where a typical application of PBC-codes would be heterogeneous catalysis, for example on metal surfaces, while gas-phase chemistry and homogeneous catalysis are usual applications of codes using localized basis sets. The largest difference however is that codes with atom-centered basis sets usually feature analytical Hessians for DFT-methods and a full set of internal coordinates for geometry optimization. This has enabled Hessian-based transition state optimization techniques such as the trust-region-image minimization,13 the PartitionedRational Function Optimization (P-RFO) approach,14 and © XXXX American Chemical Society

Received: October 27, 2017

A

DOI: 10.1021/acs.jctc.7b01070 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

In terms of the Cartesian coordinate system, the bond length bi between atoms i1 and i2 is given as

A main criticism, already early on, of constrained optimization for TS searches was that they can fail if the constraint is not chosen appropriately.18 Furthermore, the constraint has to be chosen for a specific problem and it is thus not a general approach for all kinds of transition states. Other approaches have the benefit of being, at least in principle, able to find every kind of transition state. Several methods, in particular the nudged-elastic-band (NEB) method32−36 and later the string method (SM)37 and growing string method (GSM),38−44 as well as others,45−49 directly optimize the MEP50−53 between reactant and product. This means that no knowledge besides reactant and product, e.g. no knowledge of how the reaction actually proceeds is strictly necessary. In practice, this knowledge may still be useful, and, depending on the specific method, also required, to guide these methods to convergence. The NEB method and its variants have become popular mainly for PBC-codes before, in particular SM and GSM have also been adopted in codes using localized basis functions. The dimer method54,55 has had similar success and started out being based only on gradients and Cartesian coordinates, which makes it an efficient alternative to Hessianbased optimizations methods even in codes that have analytical Hessians. Efficiency and robustness for specific problems are, besides general applicability, of course also desirable qualities for TS optimization methods. I believe that this is precisely where constrained optimization is useful although it is definitely not a general black box approach. While it is true that this requires some chemical intuition and experience, I believe that, where it works, it is the method of choice if one does have the chemical intuition and experience. For systems with PBC, constrained optimization can be extremely useful for transition states that occur in materials such as zeolites. An important aspect for zeolites is that, although the transition states one finds in acid catalysis are often not very complicated, there can be many of them. This is particularly true if one starts to explore different isomers, where the same transition state can occur at the same active site motif but in different locations of the zeolite. Here, information about the transition, as expressed in the constrained coordinate, can be reused very easily to optimize the same transition state in a different position. In this work, the implementation of constrained linear combinations of bond lengths into a geometry optimizer based on Cartesian coordinates will be discussed first. It will then be shown how, instead of a primitive scan, one can obtain the transition state in a single calculation by also optimizing the value of the constrained coordinate. Such an approach is also implemented in the molecontrol environment that is part of the Turbomole program package.5 Finally, the efficiency of this approach for a few applications to zeolite-catalyzed reactions will be demonstrated.

bi = |x i1 − x i2|

Here, x is defined as the set of all atomic Cartesian coordinates and xi1 as the subset of the three Cartesian coordinates of atom i1. The gradient of the total energy with respect to x is g = ∇E

3 = E + λ (a 0 − a)

(5)

A = ∇a

has nonzero components only for the atoms that contribute to a bond that is part of the constraint. For example, if atoms i1 and i2 form bond distance bi than the gradient with respect to the y coordinate of i1 is A xi1,y = wi(xi1, y − xi2, y)/bi

(6)

To make the Lagrangian (eq 4) stationary: ∇3 = g − λ A

(7)

=0

(8)

first of all, gradients of atoms that do not contribute to the constraint a, need to be zero. In other words, all elements of g need to be zero, for which the corresponding elements of A are zero. For the remaining elements and overall, g and A need to be linearly dependent, and then, λ is determined as λ = ATg /|A|2

(9)

In practice, one can therefore simply project the normalized gradient of the constraint a out of the energy gradient g to obtain the gradient of the Lagrangian: ⎛ AAT ⎞ ∇ 3 = ⎜1 − ⎟g |A|2 ⎠ ⎝

(10) 3

and employ this gradient in geometry optimization. In principle, this is everything that is necessary for constrained optimization with fixed combinations of bond lengths. Additionally, one may sometimes need to enforce the constraint a = a0 on a given structure x or an optimization step Δx. This may be necessary for various reasons, in the most trivial case because one wants to generate the initial Cartesian coordinates for geometry optimization with a0. There is no unique way to modify x to satisfy the constraint. In principle, all or only one involved atom can be modified in various ways to achieve a = a0. In absence of other information, the most reasonable thing to do is probably to change x in the direction of A. To first order in x, the desired change Δa with respect to the current a in Δx is therefore achieved with Δx ≈ A × Δa /|A|2

n i

(4)

where a0 is the current value of the constraint. The gradient of a

CONSTRAINED OPTIMIZATION IN CARTESIAN COORDINATES The constraint is defined as the requirement that the linear combination of a desired number n of bond lengths bi:

∑ wi × bi

(3)

The constraint can be enforced in a geometry optimization through the use of a Lagrangian:

■

a=

(2)

(11)

The Cartesian coordinates are modified according to eq 11 and the actual a is calculated. This is repeated until the desired change in Δa is converged to 10−12 Å, which typically takes less than five iterations of this process. To modify a given

(1)

equals a constant a0. The contribution of the bonds can be weighted by a factor wi. For most applications, wi is ±1. B


Article

Journal of Chemical Theory and Computation optimization step to enforce Δa = 0, one can simply project out the normalized A, as in eq 10. The derivative of E with respect a is generally not defined, since a can be changed in various ways. However, if one considers only E(a), where E(a) is assumed to be fully optimized with respect to all other coordinates, a variation δa will occur along the remaining gradient g. Since g has the same direction as A, the Cartesian variation δx for a given δa follows from eq 11. Taking the first derivative of the energy change gTδx gives

E′ =

■

d E (a) da

=g TA/|A|2

(12) (13)

SCANNING THE CONSTRAINT To determine the transition state, all one needs to do is find the local energy maximum for E(a) that corresponds to the transition state aTS. The simplest, and computationally least efficient, way to achieve this is to set up several calculation on a discrete grid of value for a. This is shown in Figure 1, where one can see how the energy varies smoothly as a function of the constrained coordinate with a maximum where the transition state is.

Figure 2. Typical 1D-optimization using (a) a constant threshold for constrained optimization and (b) an adaptive threshold. The optimization is for reaction 4 illustrated in Figure 3.

(a) Choice of minimization algorithm, not studied in this work. (b) Generation of the Cartesian coordinates of the initial structure. (c) Choice of sufficient convergence criteria. 2. Optimization of E(a) with respect to a. The choice of minimization algorithm for the constrained optimizer is of course important for the overall performance of the TS optimization. However, this is entirely independent of the method described here and one could use any optimizer, such as steepest descent, conjugate gradients, FIRE,56 or QuasiNewton methods. Here, the Quasi-Newton method with the BFGS update as implemented in the atomistic simulation environment, ASE,12 will be used. As the only modifications, the convergence criteria have been changed in case of the adaptive threshold described below. Furthermore, the optimizer is required to perform at least three steps before the optimized structure is accepted. The remaining challenges for efficient optimization listed above will now be discussed in more detail. Generation of the Cartesian Coordinates of the Initial Structure. The problem of generating a Cartesian structure for some a0 from a given structure at a was already introduced above, when discussing constrained optimization. Shifting all atoms along A generally works well. This choice can be motivated by the fact that, after a constrained optimization, A and g point in the same direction. Therefore, shifting the atoms along A also shifts them uniformly uphill in energy along g. This means that atoms that are not involved in the constraint are not modified at all. Another approach is to linearly interpolate between two Cartesian structures. This has the advantages of creating an initial structure that includes potential changes of the spectator atoms that are not part of the constraint. To prevent unrealistically large changes in these atoms, the interpolation step is only accepted when the

Figure 1. Typical scan of a constraint coordinate. The scan is for reaction 4 illustrated in Figure 3.

The transition state can then either be directly read off from a sufficiently fine grid or be redetermined by interpolating to refine aTS and performing a new optimization for aTS. An advantage of this procedure is that it can be trivially parallelized over the different, independent calculations. Alternatively, one may approach this problem as a local maximization of E in the one-dimensional space of a. This means carrying out sequential optimizations for each a to reach the TS. It is clear that, even if one performs a primitive gridscan, this will save computation time, due to better initial starting structures and overall less constrained optimizations. Figure 2a shows the outcome of a 1D-optimization of the constrained coordinate as described in more detail below. As can one see, only four steps in terms of the constrained coordinate are required to reach the transition state. This is in contrast to the scan depicted in Figure 1, that, on a much finer grid in a smaller range, yields the same level of convergence of the transition state. The following problems determine the computational efficiency of the 1D-optimization of the TS: 1. Constrained optimization of E(x) at a given value of the constraint, a0. C


Article

Journal of Chemical Theory and Computation maximum Cartesian atomic step is below the general step-limit, Δxmax, that is employed for the optimization in terms of a. Choice of Convergence Criteria. It is clear that, for a far away from aTS, one needs only relatively sloppy convergence criteria to realize that one is indeed far away from aTS and to take a reasonable step toward aTS. At aTS, however, one will want tighter convergence thresholds to obtained wellconverged transition states. Typically, the convergence of constrained optimizations is measured based on the gradient of the Lagrangian, ∇3 . When trying to find the maximum of E(a) as described above, this is not particularly useful. At a value of a that is far away from aTS, one mainly needs to determine the derivative of E(a) with respect to a, E′ with enough accuracy to obtain a better value for a. At aTS one wants to actually miminize the total gradient g. However, that a = aTS can be deduced only from a reasonably low E′. To summarize, the convergence criteria need mainly to be chosen such that E′ is accurate enough. In order to do this efficiently, the threshold for ∇3 is determined dynamically based on E′. The maximum atomic component (e.g., the maximum norm of the individual atomic force) is used to measure ∇3 : fmax (∇3) = max(|∇31| , |∇3 2|,..,|∇3 n|)

Table 1. Comparison of the Two TS Optimizations Shown in Figure 2 using 1D-Optimizations in a with an Adaptive and a Constant Threshold at Each a0a adaptive threshold

E′ × |A|2 m × fmax (A)

(14)

n(a0)

a0

f max

n(a0)

82 83 95 29 12 8

37 10 18 89 14 30 198

0.43 0.99 1.91 2.22

8 8 7 8

125 61 112 55

Δa = E′ (15)

353

Å2 eV

(19)

Step size control is applied by scaling steps back that exceed a maximum step size Δxmax. Here, Δxmax refers to the maximum atomic component of the step, e.g., the norm of the three Cartesian coordinates of the step of each atom, as in eq 14. The gradient of E(a) with respect to a, E′, is inter- or extrapolated to obtain the root, E′ = 0. If one uses a linear interpolation with two points, this corresponds to a quasiNewton step with a Hessian h obtained from the finite difference of two gradients:

(16)

Δa = −

E′ h

(20)

Thus, the initial step, eq 19, corresponds to a quasi-Newton step with a h = −1 eV/Å2. Generally, the two points with lowest E′ are used to compute the curvature, because these can be expected to provide the best estimate of the Hessian. However, care must be taken if two points are close together, because the finite difference Hessian will then become inaccurate. Especially if the convergence criterion is not tight, random numerical differences in the gradients often lead to a too large curvature for a small step. A too large curvature will result in another small step, increasing the problem. In a simple approach to prevent this, an average curvature is computed by linear regression of all E′(a) vs a. A two-point finite difference curvatures is restricted to deviate at most by a factor of 10 from the average Hessian. If the root of E′ is in the interval of two previous calculations, a1 and a2, as determined from an opposite sign of E′(a1) and E′(a2), steps should generally be within this interval. Inaccurate gradients may, in the worst case, have a wrong sign, thus putting wrong limits on the position of the root. Eventually, such a wrong interval for the root will result in an optimization step that violates this interval, going toward the true position of the root. However, a step violating an interval from previous steps may also be due to an inaccurate Hessian leading to a too large step. If the interval is correct, one may, instead of a quasi-

(17)

with gmax = 0.01 eV/Å. The optimization with an adaptive threshold will therefore stop when either eq 16 or eq 17 is fulfilled. To reliably achieve this convergence criterion in eq 17 if no adaptive threshold is used, the constant threshold for 3 is set to fmax (∇3) ≤ 0.8 × gmax

f max

dimensional, the direction follows from the sign of the gradient E′ and only the step size needs to be determined. The initial step as well the step in the case of a wrong curvature of the energy is determined from the gradient as

Additionally, fmax (∇3) is always to required to be ≤0.1 eV/Å. The final convergence criterion for the transition state employed in this work is

fmax (g) ≤ gmax

a0 0.43 0.99 1.90 2.31 2.25 2.23

Columns 2−4 show the optimization with an adaptive threshold, and columns 5−7 show the optimization with a constant threshold. The column labeled a0 shows the value of the constrained coordinate in Å at the respective step. Column f max lists fmax (∇3) in meV/Å, and the column labeled n(a0) refers to the number of steps used in the constrained optization at a0.

where m is the number of atoms that constitute the bond distances present in the constraint according to eqs 1 and 2. It is not possible to set a threshold for fmax (∇3) such that it would guarantee a certain accuracy of E′. However, following eq 15, the convergence of fmax (∇3) is controlled with a relative threshold f thr, which is typically around 0.4 to 0.8: fmax (∇3) ≤ fthr ×

step 0 1 2 3 4 5 total a

It follows from eqs 13 and 14 that E′ ≤ m × fmax (g) × fmax (A)/|A|2

constant threshold

(18)

Optimizations using a constant and an adaptive threshold are compared in Figure 2 and Table 1. Due to the more accurate gradient, E′, the optimization with a constant threshold requires less optimization steps in terms of the constrained coordinate, a (4 vs 6). However, this comes at the cost of more energy/gradient evaluations (353 vs 198). The lower level of convergence far away from the TS is noticeable from the gradients listed in Table 1. It is also visible from the energy curve in Figure 2. Optimization of E(a) with Respect to a. The root of E′ is determined using a combination of interval search and quasiNewton algorithm. Since the optimization is only oneD


Article

Journal of Chemical Theory and Computation Newton step, proceed with an interval-search step, e.g., compute a point right in between the values a1 and a2 that bracket the root and then resume the usual quasi-Newton type optimization. If the interval is not correct, one can increase the thresholds used for constrained optimization to obtain more accurate gradients. Since it is, without further calculation, impossible to tell whether the interval is correct or not, the following strategy is adopted. When a quasi-Newton step violates an interval, interval steps can be performed at most three times in total which should be sufficient to find the root, if the interval is correct. When a step violates the interval a fourth time, the entire transition state search is automatically restarted from the current structure with new convergence criteria. If an adaptive convergence criterion has been used, the new search is performed with a constant convergence criterion using the default according to eq 18. If the previous optimization has already been performed with a constant convergence criterion (>0.1 × gmax) the criterion is reduced to (0.1 × gmax). Otherwise the optimization is terminated unsuccessfully. Generally, restarting with tighter criteria from the current structure does not lead to a too dramatic increase in computation time since violation of an interval will usually occur very close to the TS. In the performed test calucations, this restart was only very seldom necessary. However, it generally combines the efficiency of using an adaptive threshold with the increased robustness of a constant threshold that is sometimes necessary. A quadratic interpolation of three gradients was also tested. This is equivalent to including third-order (cubic) anharmonic contributions in terms of E(a). The use of a cubic step is regulated with a threshold that accepts the cubic step only if it does not deviate from the quasi-Newton step significantly. The ratio of the cubic step to the uasi-Newton step must not exceed this threshold, which in the test calculation was set to 4.0.

■

BENCHMARK CALCULATIONS Systems Used for Testing. Several test calculations were run that serve to benchmark the efficiency of the method and allow the parameters employed in the optimization to be optimized. All tests used for transition state optimization were taken from real problems that were actually optimized using the described method in a recent investigation of the initiation of the MTO process for the zeolite H-SSZ-13 (CHA).57 The reactions, a sketch of the transitions state as well as the constrained coordinate are depicted in Figure 3. The classic example of constrained optimization for TS optimization is an SN2-type reaction, where a group is transferred between two molecules. For zeolites, this type of reaction is important for a wide range of reactions such as protonations and alkylations. Since all of these cases are relatively similar and since the described method works well, only three methylation reactions are considered. Here the methylation of methanol, water, and propylene is considered as important steps in the MTO process. Additional methylation reactions that have been studied successfully in ref 57 are methylations of various olefins, CO, and substituted ketenes. Furthermore, many alkylation reactions involving formaldehyde, trioxane, dimethoxymethane, and other acetals have been investigated with this approach.58 The other reactions depicted in Figure 3 involve more bonds and therefore highlight how more complicated reactions can be tackled with this method. Reactions 4 and 5 involve the dehydrogenation of MeOH and formaldehyde. In these reactions, the acid site protonates a C−H bond yielding H2, while another O−H or C−H bond is

Figure 3. Transition states studied in this work. Panel a shows the involved reactions and panel b shows an illustration of the transition states studied in this work. Involved formed and broken bonds are highlighted in red. The constrained coordinates are also shown.

simultaneously deprotonated. In reaction 6, methyl formate is decomposed to methanol and CO. All of the transition states so far can be dealt with by simply constructing the constrained coordinate using all of the bonds that are formed/broken with weights of ±1. The next TS is special in the sense that this straightforward assignment does not work. Reaction 7 is part of the methane−formaldehyde mechanism that has been proposed to be responsible for the initiation mechanism of the MTO process in the absence of impurities. Here, methane is deprotonated and the CH3−-moiety performs a nucleophilic attack at formaldehyde, which is then protonated at the oxygen to form ethanol. It turns out that formaldehyde is relatively easily protonated in the initial structure and is then only loosely attached to the acid site. If one includes the O-protonation of formaldehyde in the constraint (with the usual ±1 weights), E


Article

Journal of Chemical Theory and Computation Table 2. Results of the TS Optimizations Using Different Parameters and the Bad Starting Guessa f thr

0.4

0.6

0.6

0.6

0.8

-

-

Δxmax

0.2

0.2

0.2

0.2

0.2

0.2

-

-

cubic

-

-

4.0

-

-

-

-

-

interp.

-

-

-

T

-

-

-

-

dimer

-

-

-

-

-

-

0.01

0.05

1 2 3 4 5 6 7 8 9 total mean median

191 160 283 217 188 364 224 223 647 2496 277 223

165 218 288 201 186 281 251 210 551 2351 261 218

166 174 287 255 217 305 246 201 598 2448 272 246

166 169 288 222 178 355 192 206 513 2289 254 206

145 169 278 196 215 283 299 257 703 2544 283 257

251 192 356 353 230 *618 339 358 572 3269 363 353

488 440 1261 *1471 702 1694 -

465 397 1269 699 -

-

a

The top part of the table describes the method used for optimization, and the bottom part lists the number of energy/gradient evaluations for the reactions under study as well as the sum of these. Columns 2−7 list different variations of 1D-optimization using constraints and columns 8 and 9 list calculations with the dimer method using two different offsets between the dimer (Å), given as the number in the row labeled dimer. f thr denotes the use of an adaptive threshold according to eq 16, Δxmax is the maximum atomic Cartesian step size, cubic is the threshold for the use of a cubic optimization step, and interp. is whether or not linear interpolation between Cartesian structures is used for the generation of initial structures. A “-” signifies that a feature has not been used at all. The result of calculations that did not converge, using the described settings, is given as “-”. Calculations that converged to a different conformer of the same TS are labeled with an “*”.

all opimitization have been repeated three times, and the average value is used to evaluate the performance. The variation is usually less than 10%. All optimizations were carried out using a modified version of the ASE12 (version 3.10.0), where the optimizer was implemented. As already mentioned the constrained optimization for a given value a0, an energy minimization, will be carried out with the BFGS optimizer of the ASE. It is important to realize that the number of gradient evaluations depends both on the underlying optimizer and its implementation (BFGS, ASE) as well as the specific system (zeolite, CHA). Using a different optimizer or a different (larger) zeolite could systematically increase or decrease the performance. The transition states optimized in this work have typically around 120 atoms. As a reference point, a typical full minimization usually takes around 100 steps, in some cases more. For example, optimization of an isomer of the acid site (Al−OH) removing the proton from one oxygen and placing it at the other oxygen atom took 70 iterations. Optimization of a minimum with a SMS (Al-OMe) required 183 steps using a preoptimized CH3 group. Optimizing the hydrogen bound adsorbates of methanol and methyl formate at the acid site (Al−OH) took 150 to 250 steps depending on initial guess and isomer. As will be shown below, a TS optimization with the described method typically takes around 200 to 300 steps, e.g., it is about two to three times slower than a minimization. A transition state that involves a hydrogen bound adsorbate required more than 500 steps, thus being slower as observed for the related minimizations mentioned above. I expect that this ratio will roughly stay the same when using a different optimizer or zeolite. In particular, if one were to use a more efficient optimizer, I expect that both minimizations and TS optimization would benefit in the same way. In different systems with more/less degrees of freedom both minimization and TS optimization can be expected to require more/less steps.

variation of the value of the constraint leads to initial protonation followed by dissociation of the protonated formaldehyde from the active site, while the desired reaction does not occur. If one omits the O-protonation of formaldehyde in the constraint and scans the remaining constraint, formaldehyde is simply initially protonated and then reacts in the desired way. This highlights the challenges with the definition of the constrained coordinate, even in case where it works well, once properly defined. Reaction 8 is equivalent to reaction 4 in the sense that methanol is oxidized to formaldehyde; however, here methane is formed instead of hydrogen. Reaction 9 is similar to reaction 1, only here, methylation of methanol occurs in the direct pathway,59 by an adsorbed methanol, rather than the indirect one, by a previously generated SMS. Benchmark Procedure. The performance will be measured in terms of the number of gradient evaluations to reach the convergence criterion (all atomic forces below gmax = 0.01 eV/Å). Energy and gradients are always calculated both, e.g., no attempt is made to save computation time by evaluating only the energy and not the gradient in certain steps. Furthermore, another approach, such as machine-learning,60 to generally reduce the cost of geometry optimizations, could be used but has not been explored in this work. Energies and gradients are calculated with density functional theory (PBE-D361,62) and the projector-augmented-wave (PAW) method using the VASP program package in version 5.4 as well as the standard VASPPAWs.63,64 The Brillouin-zone was sampled at the Γ-point and Gaussian-smearing with a width of 0.1 eV was used. An energycutoff of 400 eV was used for the expansion of the wave function in the plane wave basis set. The initial wave function is randomized. The convergence criterion used for optimization of the constraint is low far away from the TS when an adaptive threshold is employed. Together with the initial randomization of the wave function, this leads to some variations in terms of the steps required to reach convergence. As this is a real effect, F


Article


Table 3. Results of the TS Optimizations for the Automated Relaxed Potential Energy Surface Scan, the Dimer Method, and the CI-NEB Method Using a Good Starting Guessa scan guess 1 2 3 4 5 6 7 8 9

+10 73 60 75 89 140 97 153 99 126

−10 68 57 80 119 102 73 153 102 129

±10 33 25 24 53 41 61 75 28 25

dimer +20 106 118 95 172 185 212 188 146 189

−20 98 81 112 198 187 110 160 105 182

±20 54 49 37 49 51 59 86 65 74

+10 110 120 128 177 261 294 83 -

−10 114 219 204 499 105 144

±10 67 79 63 71 87 59 75 51

1-CI-NEB

+20 153 272 369 249 297 230 348

−20 137 284 207 398 -

±20 99 111 111 115 159 107 127 143 99

±10 13 14 18 22 23 28 18 23 16

±20 774 459 29 46 46 151 132 721

3-CI-NEB ±10 42 45 63 135 105 165 75 81 66

±20 156 261 123 144 192 210 318

The row labeled “guess” specifies the information used, e.g., a structure removed from the transition state by 10 optimization steps in either direction (+10, −10) or both (±10) and analogously for 20 steps. The direction of displacement (+ or -) is ambigous and is just used to label these structures. The number of optimized images for the CI-NEB is specified as the leading integer. The result of calculations that did not converge, using the described settings, is given as “-”. a

CH3 moiety away from the active site and constructing a planar CH3-group would clearly be a better starting point. The initial structures are neither a minimum nor very close to the transition state, as they have on average four to five sizable harmonic imaginary frequencies. Table 2 lists the number of energy/gradient evaluations needed to reach convergence for the nine studied reactions. The main conclusion is that calculation with a fixed convergence threshold is generally about 50% less efficient than using the adaptive thresholds. The influence of other parameters, if one uses the adaptive threshold, are relatively small. Values of f thr > 1 often caused serious convergence problems, as one might expect. Values of f thr < 0.4 generally lead to increased computational cost due to unnecessarily accurate constrained optimizations distant from the TS. For values 0.6 ≤ f thr ≤ 0.8, only small variations in the performance were found, since the cost and benefits of an accurate gradient E′ are to some extent balanced. This is because a less tight gradient criterion will save time in terms of the number of steps per constrained optimization but may lead to more optimizations steps in terms of the constrained coordinate due to a lower accuracy of the gradient with respect to the constraint. Using a cubic rather than a quadratic optimization scheme did also not lead to significant improvements. The same is true for the linear interpolation procedure, which does not improve or worsen the results. When an adaptive threshold is used, typically 200−300 iterations are required to obtain the converged transition state using the starting guess described above. This means that TS calculation is 2 to 3 times more expensive then minimization. The dimer method calculations did not converge in all cases, which is largely due to the bad starting guess that does not have the correct Hessian eigenvalue spectrum. Where the dimer method converged, it required about 3−5 times more gradient calls than TS optimization based on the constraints, e.g., 4 to >10 times the cost of a minimization. The distances between the dimers (0.01 Å) and the finite difference steps for dimer rotation (dθ = 10−4 rad and dR = 10−3 Å) were kept at the default values as in the original publication of the dimer method, ref 54. The dimer distance corresponds to the offset in a finite difference Hessian calculation. Given the strong curvature of the transition states and a SCF-convergence criterion of 10−9 eV, the choice of separation is expected to be

Of the methods one could compare to, the dimer method is the most obvious. This is because it also requires only a single reasonable starting guess. Using the implementation of the dimer method in the ASE, a calculation can be initialized with an initial search direction, or the dimer method will determine a search direction based on a random initial guess. Using a random guess rarely worked, e.g., the corresponding dimer method calculations did usually not converge. This can be explained by the relatively bad starting guess employed here. The initial Hessian (which is not used in any of the optimizations) has generally not the right eigenvalue spectrum. Supplying the dimer method with an initial search direction (A) that is equivalent to the constraint often led to convergence and this is the approach that has been used in all cases discussed in this work. Importantly, if one initializes the dimer method with the specific search direction, it takes the same user time to set up as a constrained optimization. In this case, the main advantage of the dimer method is that the search direction is not limited to the initial direction. Therefore, in cases where a constrained optimization fails, the equivalent dimer method calculation using the same information may still succeed. Results Using a Bad Starting Guess. The performance of an optimization procedure, and in particular TS optimization, depends sensitively on the starting guess. To illustrate this, all transition state searches were carried out using both a bad and a good starting guess. A bad starting guess was constructed using a “copy and paste” approach that is intended to be unbiased by knowledge about the optimized TS and not significantly better than one can expect in a real application. In fact, the construction procedure gives initial guesses that are generally worse than they could be when employing some chemical intuition. All TS involve the reaction of one or two molecules with the acid site of a zeolite (a bridging surface hydroxyl group) or a SMS (a bridging methoxy group). The initial guess was constructed by manually placing the separately preoptimized molecules in the preoptimized zeolite, without altering their individual structures. The placement follows the same chemical intuition that is also required in determining the involved bonds. For example, in the reaction where methanol is methylated, the oxygen of methanol is placed in front of the methoxy group with a C−O distance of 2.0 Å, forming roughly one line from the methanoloxygen, the CH3-carbon, and the zeolite-oxygen. However, the CH3 moiety of the SMS is not modified, although moving the G


Article


naturally provide well-defined end points for the CI-NEB method.

appropriate. An additionally tested dimer separation of 0.05 Å did not lead to a large difference in performance. Results Using a Good Starting Guess. A good starting guess was obtained by distorting the optimized transition state in both directions along the transition mode and optimizing the obtained structure, thus roughly following the reaction path in both directions. Guesses of different quality can now be obtained by controlling the number of optimization steps. If one performs a full optimization of reactant and product, the Cartesian distance to the transition state is rather large. This is because reactant and product typically contain loosely bound molecules that are rotated and translated significantly. While these rotations and translation are chemically insignificant, they complicate any transition optimization procedure, especially one that employs Cartesian coordinates. The Cartesian distance between fully optimized reactants and products is typically around 30 Å (always >15 Å), while it is ∼4 Å after 10 optimization steps and ∼7 Å after 20 optimization steps. In the benchmark calculations, an initial guess based on both 10 and 20 optimization steps in both directions was used. This is likely a better guess than one could usually construct, without knowing the transition state beforehand or carrying out preoptimizations of some kind. The calculations were now only performed for the best setting for the automated relaxed potential energy surface scan (f thr = 0.6, Δxmax = 0.2 Å) and for a dimer separation of 0.01 Å. The results are shown in Table 3. Both methods show good performance and converge in around 100 or sometimes less iterations. Since now both end points are available, one can directly compare to the climbing-image nudged elastic band method (CI-NEB33,34), which employs both structures (labeled ±10). As only short paths are optimized here, a small number of optimized images (in addition to the two end points) is used: Three and one image, e.g., a total of five and three images. Within a single calculation, a NEB was first run to a convergence criterion of 0.1 eV/Å followed by a CI-NEB optimization to the full convergence criterion of 0.01 eV/Å, using the implementation in the ASE12 (version 3.10.0) with default parameters. The CI-NEB approach works very well, when three optimized images are used and requires less than 30 iterations when only a single image is used. This is partially due to the very short reaction path, where the initial linear interpolation already provides a very good guess for the transition state. If one also starts the automated relaxed potential energy surface scan and the dimer method calculation from a central interpolation between the two end points (labeled ±10 in Table 3) the results are improved, too. For the dimer method, the initial direction has in this case been obtained as the connecting vector between initial and final state, which gave an additional, slight improvement over the usual initialization using the predefined bond lengths (latter results not shown). Additional CI-NEB calculations using a longer reaction path (±20, ∼7 Å) are also shown in Table 3 and require around a few hundred steps when three optimized images are used. Using fully optimized end points significantly increases the computational cost. Overall, all three tested methods work well when provided with good starting guesses. Which method one prefers will also depend on which information one finds easiest to provide, reasonably close end points or a TS-guess with an initial mode or broken and formed bonds. This also varies with the systems which one studies. For example, reactions on metal surfaces often lead to tightly bound reactants and products which

■

SUMMARY AND CONCLUSION Geometry optimization based on Cartesian coordinates was combined with constrained linear combinations of bond lengths. Transition states can be obtained through multiple constrained optimizations by scanning along the constraint. This can be accomplished in a single calculation by finding the root of the gradient with respect to the constraint. Such an optimization is most efficient when an adaptive threshold is used for the constrained optimization. The threshold depends on the gradient with respect to the constraint; that is, the threshold will become tighter when the optimization approaches the transition state. Other details of the optimization algorithm were found to have only little influence on the performance. Generally, the optimizer used for constrained optimization will have a large impact on the overall performance. This optimizer can be changed trivially; however, this has not been explored in this work, and it is expected that the effect is similar to that for usual energy minimizations. This approach for TS optimization has already been used for zeolite-catalyzed reaction, and a few test cases have been reevaluated in this work. Transition states are obtained reliably from a readily constructed and not particularly good initial guess. A TS optimization is about two to three times more costly than a typical minimization. Scanning coordinates is not a black box approach to transition state optimization, and it does not work in all cases. Setting it up correctly and anticipating when it works and when it does not requires experience. However, for the type of reactions where it works, it is very robust and not very sensitive to the starting guess. For reactions that are well-defined in terms of broken and formed bonds, for example in zeolite catalysis, this approach is applicable to many reactions. When the same type of transition state needs to be reoptimized for a different isomer, as is often the case in zeolite catalysis, information about the structure of the TS can be efficiently reused through the value of the constrained coordinate. Overall, I believe that constrained optimization using linear combinations of bond lengths is a useful tool that is often complementary to established approaches (NEB, dimer, SM, and GSM) when optimizing transition states in periodic structures.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b01070. Cartesian coordinates of input structures and optimized transition states (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Philipp N. Plessow: 0000-0001-9913-4049 Funding

The author acknowledges support by the state of BadenWürttemberg through bwHPC (bwunicluster and JUSTUS, RV H


Article


(22) Bell, S.; Crighton, J. S. Locating transition states. J. Chem. Phys. 1984, 80, 2464−2475. (23) Quapp, W.; Hirsch, M.; Imig, O.; Heidrich, D. Searching for saddle points of potential energy surfaces by following a reduced gradient. J. Comput. Chem. 1998, 19, 1087−1100. (24) Quapp, W.; Hirsch, M.; Heidrich, D. Bifurcation of reaction pathways: the set of valley ridge inflection points of a simple threedimensional potential energy surface. Theor. Chem. Acc. 1998, 100, 285−299. (25) Munter, T. R.; Bligaard, T.; Christensen, C. H.; Nørskov, J. K. BEP relations for N2 dissociation over stepped transition metal and alloy surfaces. Phys. Chem. Chem. Phys. 2008, 10, 5202−6. (26) Liu, X.; Xiao, J.; Peng, H.; Hong, X.; Chan, K.; Norskov, J. K. Understanding trends in electrochemical carbon dioxide reduction rates. Nat. Commun. 2017, 8, 15438. (27) Vojvodic, A.; Calle-Vallejo, F.; Guo, W.; Wang, S.; Toftelund, A.; Studt, F.; Martinez, J. I.; Shen, J.; Man, I. C.; Rossmeisl, J.; et al. On the Behavior of Brønsted-Evans-Polanyi relations for Transition Metal Oxides. J. Chem. Phys. 2011, 134, 244509. (28) Kudin, K. N.; Scuseria, G. E.; Schlegel, H. B. A redundant internal coordinate algorithm for optimization of periodic systems. J. Chem. Phys. 2001, 114 (7), 2919−2923. (29) Bucko, T.; Hafner, J.; Angyan, J. G. Geometry optimization of periodic systems using internal coordinates. J. Chem. Phys. 2005, 122 (12), 124508. (30) Bucko, T.; Benco, L.; Dubay, O.; Dellago, C.; Hafner, J. Mechanism of alkane dehydrogenation catalyzed by acidic zeolites: Ab initio transition path sampling. J. Chem. Phys. 2009, 131 (21), 214508. (31) Jafari, M.; Zimmerman, P. M. Reliable and efficient reaction path and transition state finding for surface reactions with the growing string method. J. Comput. Chem. 2017, 38, 645−658. (32) Mills, G.; Jonsson, H. Quantum and thermal effects in H2 dissociative adsorption: evaluation of free energy barriers in multidimensional quantum systems. Phys. Rev. Lett. 1994, 72, 1124−1127. (33) Henkelman, G.; Uberuaga, B. P.; Jonsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys. 2000, 113, 9901−9904. (34) Henkelman, G.; Jonsson, H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. J. Chem. Phys. 2000, 113, 9978−9985. (35) Trygubenko, S.; Wales, D. A doubly nudged elastic band method for finding transition states. J. Chem. Phys. 2004, 120, 2082− 2094. (36) Bohner, M. U.; Meisner, J.; Kastner, J. A QuadraticallyConverging Nudged Elastic Band Optimizer. J. Chem. Theory Comput. 2013, 9, 3498−504. (37) Weinan, E.; Ren, W.; Vanden-Eijnden, E. String method for the study of rare events. Phys. Rev. B: Condens. Matter Mater. Phys. 2002, 66, 052301. (38) Peters, B.; Heyden, A.; Bell, A.; Chakraborty, A. A growing string method for determining transition states: Comparison to the nudged elastic band and string methods. J. Chem. Phys. 2004, 120, 7877−7886. (39) Goodrow, A.; Bell, A. T.; Head-Gordon, M. Transition statefinding strategies for use with the growing string method. J. Chem. Phys. 2009, 130, 244108. (40) Goodrow, A.; Bell, A. T.; Head-Gordon, M. Development and application of a hybrid method involving interpolation and ab initio calculations for the determination of transition states. J. Chem. Phys. 2008, 129, 174109. (41) Behn, A.; Zimmerman, P. M.; Bell, A. T.; Head-Gordon, M. Incorporating linear synchronous transit interpolation into the growing string method: Algorithm and applications. J. Chem. Theory Comput. 2011, 7, 4019−4025. (42) Goodrow, A.; Bell, A. T.; Head-Gordon, M. A strategy for obtaining a more accurate transition state estimate using the growing string method. Chem. Phys. Lett. 2010, 484, 392−398.

bw16G001 and bw17D011). Financial support from the Helmholtz Association is also gratefully acknowledged. Notes

The author declares no competing financial interest. The source code will be made publicly available at https:// github.com/pplessow/scan-ts.

■

REFERENCES

(1) Sheppard, D.; Terrell, R.; Henkelman, G. Optimization methods for finding minimum energy paths. J. Chem. Phys. 2008, 128, 134106. (2) Schlegel, H. B. Exploring Potential Energy Surfaces for Chemical Reactions: An Overview of Some Practical Methods. J. Comput. Chem. 2003, 24, 1514−1527. (3) Schlegel, H. B. Geometry optimization. Wiley Interdisciplinary Reviews: Computational Molecular Science 2011, 1, 790−809. (4) Frisch, M. J.; et al. Gaussian 16, revision A.03; Gaussian Inc.: Wallingford, CT, 2016. (5) TURBOMOLE V7.2 2017, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989−2007, TURBOMOLE GmbH, since 2007; available from http://www. turbomole.com. (6) Shao, Y.; Gan, Z.; Epifanovsky, E.; Gilbert, A. T. B.; Wormit, M.; Kussmann, J.; Lange, A. W.; Behn, A.; Deng, J.; Feng, X.; et al. Advances in molecular quantum chemistry contained in the Q-Chem 4 program package. Mol. Phys. 2015, 113, 184−215. (7) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347−1363. (8) Kresse, G.; Hafner, J. Ab initio molecular dynamics for liquid metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 47, 558−561. (9) Kresse, G.; Hafner, J. Ab initio molecular-dynamics simulation of the liquid-metal-amorphous-semiconductor transition in germanium. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 49, 14251−14269. (10) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter 2009, 21, 395502. (11) Mortensen, J. J.; Hansen, L. B.; Jacobsen, K. W. Real-space grid implementation of the projector augmented wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 71, 035109. (12) Larsen, A. H.; Mortensen, J. J.; Blomqvist, J.; Castelli, I. E.; Christensen, R.; Dułak, M.; Friis, J.; Groves, M. N.; Hammer, B.; Hargus, C.; et al. The atomic simulation environment-a Python library for working with atoms. J. Phys.: Condens. Matter 2017, 29, 273002. (13) Helgaker, T. Transition-state optimizations by trust-region image minimization. Chem. Phys. Lett. 1991, 182, 503−510. (14) Baker, J. An Algorithm for the Location of Transition States. J. Comput. Chem. 1986, 7, 385−395. (15) Ayala, P.; Schlegel, H. A combined method for determining reaction paths, minima, and transition state geometries. J. Chem. Phys. 1997, 107, 375−384. (16) Peng, C.; Ayala, P. Y.; Schlegel, H. B.; Frisch, M. J. Using redundant internal coordinates to optimize equilibrium geometries and transition states. J. Comput. Chem. 1996, 17, 49−56. (17) Fukui, K. A Formulation of the Reaction Coordinate. J. Phys. Chem. 1970, 74, 4161. (18) Williams, I. H.; Maggiora, G. M. Use and Abuse of the Distinguished-Coordinate Method for Transition-State Structure Searching. J. Mol. Struct.: THEOCHEM 1982, 89, 365−378. (19) Rothman, M. J.; Lohr, L. L. Analysis of an energy minimization method for locating transition states on potential energy hypersurfaces. Chem. Phys. Lett. 1980, 70, 405−409. (20) Scharfenberg, P. An Improved Method for the Evaluation of Transition-States. J. Comput. Chem. 1982, 3, 277−282. (21) Scharfenberg, P. Theoretical-Analysis of Constrained Minimum Energy Paths. Chem. Phys. Lett. 1981, 79, 115−117. I


Article

Journal of Chemical Theory and Computation (43) Zimmerman, P. M. Growing string method with interpolation and optimization in internal coordinates: method and examples. J. Chem. Phys. 2013, 138, 184102. (44) Zimmerman, P. M. Automated discovery of chemically reasonable elementary reaction steps. J. Comput. Chem. 2013, 34, 1385−92. (45) Burger, S.; Yang, W. Quadratic string method for determining the minimum-energy path based on multiobjective optimization. J. Chem. Phys. 2006, 124, 054109. (46) Behn, A.; Zimmerman, P. M.; Bell, A. T.; Head-Gordon, M. J. Chem. Phys. 2011, 135, 224108. (47) Plessow, P. Reaction Path Optimization without NEB Springs or Interpolation Algorithms. J. Chem. Theory Comput. 2013, 9, 1305. (48) Halgren, T.; Lipscomb, W. The synchronous-transit method for determining reaction pathways and locating molecular transition states. Chem. Phys. Lett. 1977, 49, 225−232. (49) Elber, R.; Karplus, M. A method for determining reaction paths in large molecules: Application to myoglobin. Chem. Phys. Lett. 1987, 139, 375−380. (50) Gonzalez, C.; Schlegel, H. An improved algorithm for reaction path following. J. Chem. Phys. 1989, 90, 2154−2161. (51) Gonzalez, C.; Schlegel, H. Reaction path following in massweighted internal coordinates. J. Phys. Chem. 1990, 94, 5523−5527. (52) Birkholz, A. B.; Schlegel, H. B. Path optimization by a variational reaction coordinate method. I. Development of formalism and algorithms. J. Chem. Phys. 2015, 143, 244101. (53) Birkholz, A. B.; Schlegel, H. B. Path optimization by a variational reaction coordinate method. II. Improved computational efficiency through internal coordinates and surface interpolation. J. Chem. Phys. 2016, 144, 184101. (54) Henkelman, G.; Jonsson, H. A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives. J. Chem. Phys. 1999, 111, 7010−7022. (55) Kastner, J.; Sherwood, P. Superlinearly converging dimer method for transition state search. J. Chem. Phys. 2008, 128, 014106. (56) Bitzek, E.; Koskinen, P.; Gahler, F.; Moseler, M.; Gumbsch, P. Structural relaxation made simple. Phys. Rev. Lett. 2006, 97, 170201. (57) Plessow, P. N.; Studt, F. Unraveling the Mechanism of the Initiation Reaction of the Methanol-to-Olefins Process using Ab Initio and DFT calculations. ACS Catal. 2017, 7, 7987. (58) Goncalves, T. J.; Arnold, U.; Plessow, P. N.; Studt, F. Theoretical Investigation of the Acid Catalyzed Formation of Oxymethylene Dimethyl Ethers from Trioxane and Dimethoxymethane. ACS Catal. 2017, 7, 3615−3621. (59) Ghorbanpour, A.; Rimer, J. D.; Grabow, L. C. Computational Assessment of the Dominant Factors Governing the Mechanism of Methanol Dehydration over H-ZSM-5 with Heterogeneous Aluminum Distribution. ACS Catal. 2016, 6, 2287−2298. (60) Peterson, A. A. Acceleration of saddle-point searches with machine learning. J. Chem. Phys. 2016, 145, 074106. (61) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1997, 78, 1396−1396. (62) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (63) Kresse, G.; Furthmuller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169−11186. (64) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 1758−1775.

J


Efficient Transition State Optimization of Periodic Structures through

Recommend Documents