Reaction Path Optimization

5.8 Reaction Path Optimization

5.8.1 Background and Program structure

The goal of self-consistent optimization of the reaction path (RP) is usually to obtain an initial guess for Transition State Search or an approximation to the barrier. Methods that use reactant and product structure to compute the RP are often referred to as ’double-ended’ methods or, if the the RP is discretized, ’chain-of-states’ methods. [50–52]

The RP connects reactant and product, its highest point being the transition state. It is a steepest descent path, which means that its tangent t is always parallel to the gradient g. The RP is in actual calculations discretized by a finite number of structures n. The tangents are parallel to the gradients for all structures i = 1,..,n. Assuming normalized tangents t_i^Tt_i = 1, this can be written as:

0 = (1 - t_it_i^T)g_i

(5.6)

Several approximations for the tangents t_i exist [52,53], usually using finite difference schemes. The most common methods, the Nudged Elastic Band (NEB) [52,53] and String Method (SM) [54] prevent the structures from ’sliding’ down the reaction path towards products and reactants with additional springs or interpolation/redistribution algorithms. The method used here achieves equal spacing via constrained optimization assuming a quadratic potential. [55] An initial path is provided using a slight variation of the Linear Synchronous Transit [50].

The structure of the optimization is the same as in other TURBOMOLE structure optimizations: As the jobex script drives optimizations by calling statpt/relax as well the SCF and gradient modules.The woelfling-job script drives optimizations by calling woelfling as well as the SCF and gradient modules. The woelfling-job scripts reads the current path from file path.xyz which is the output of the woelfling program. woelfling-job then creates folders to run the calculations of each structure in, gathers coordinates and gradients, and then calls woelfling again.

The aim of RP optimization is usually not to optimize the RP to some accuracy, but to obtain an initial guess for a TS optimization. It is in general not possible to find a convergence criterion (and a corresponding threshold) that guarantees a good initial guess. The maximum number of iterations and the convergence threshold are therefore relatively high and tight. One can extract a TS guess also during the course of the optimization. If the TS search is successful (or not) you can stop (or restart) the RP optimization. Apart from simply killing the program you can add a ’stop’ file in the (scratch) directory, in which the script runs. It will then terminate at the end of the current cycle and can easily be restarted.

5.8.2 Input Structure

Options can be modified using keywords in the $woelfling data-group. The most important options are:

$woelfling
ninter           14
riter 0
ncoord            2
align             0
maxit          40
dlst   3.00000000000000
thr  1.000000000000000E-004
method q

The values above are the default values. If $woelfling is missing, it will be added during the first woelfling run and default values will be set. Most importantly, ncoord is the number of input structures provided, ninter is the number of structures to be used for discretization of the path and maxit is the number of cycles to run. If align 0, structures will be rotated/translated to minimize the cartesian distance, for align 1 strutures will be used as provided. Using method qg instead of method q, a reaction path will be grown as in the growing string method. To start a RP optimization you need to provide at least a reactant and a product structure ncoord ≥ 2. You may provide more structures if you have a guess for the reaction path. The input structures will be used to compute an initial guess with ninter structures that is then optimized. Reactant and product structure will stay fixed throughout the optimization. All structures have to have the same ordering of atoms!

The input structures are provided in a file coords, which contains merged coord files. All ncoord structures are given in the right order in the typical TURBOMOLE coord format.

5.8.3 How it works

Minimum Input/Quick and Dirty

1.: Make a usual TURBOMOLE input using the coord file of either reactant of product structure.
2.: Join the coord-files of reactant and product in a file coords.
3.: Run woelfling-job
4.: Check the output and the path (path.xyz) to extract a TS guess.

It is usually a good idea to check the initial path before starting the calculation. Once you have prepared the input, simply run woelfling directly and check path.xyz. If it looks reasonable, just run the woelfling-job script.

Unsuccessful Optimization If there is no reasonable TS guess or a frequency calculation does not give the correct number of imaginary frequencies, you can:

1.

Check if have used the best structure as TS guess, maybe the structure with the highest energy is not the best.

2.

Check if the RP has a reasonable amount of structures (if they are far apart, it is unlikely that a structure is close to the TS)

3.

Check if the RP is reasonably converged (mean of rms(g_i^⊥) in output <1.0d-3; path is continuous in terms of energy and structure)

(a): If it is not yet converged, converge it.
(b): If it is not going to converge, provide useful structures for the initial guess or maybe use more structures for the path.

Restart

1.: If you have stopped the calculation adding a ’stop’ file, you can just run woelfling-job again.
2.: If you have run the maximum number of cycles, just increase maxit and run woelfling-job again. It will then run from the old maxit to the new maxit.

If the calculation crashed ’on its own’ it is likely that the SCF failed to converge - improve the corresponding options. If the optimization crashed in the middle of the run, it is most likely that it crashed during a SCF or gradient step, since basically all cpu time is spent there. In that case, remove at least the folder where the SCF and gradient program had been running when the program crashed. The files necessary for woelfling should be intact and you can simply restart it.

Restart - more details If the files are not intact, one can still use the optimized coordinates in one way or the other. There are basically three phases which are explicitely indicated by the number of iterations riter =:

0: woelfling reads control, coords to generate an initial guess path.xyz
1: woelfling reads control, gradients to compute an optimized path.xyz. Hessian are initialized and written to hessians-new
>1: woelfling reads control, gradients, oldgradients and hessians to update hessians and compute an optimized path.xyz. Hessian are updated and written to ’hessians-new’.

Therefore repeated execution of woelfling will yield the same output (unless new energies/gradients are computed). The woelfling-job script takes care of the file-handling and should enable restart at any time. If files are damaged it will be hard to gather gradients and corresponding hessian- and oldgradient-information. If you have an intact gradients or oldgradients-file (necessary condition: number of lines=(3 + 2 × natoms)×number of structures), name it gradients, set riter in the control-file to 1 and restart. If those file are not intact, you can extract whatever structure information you have obtained and use it to provide a better initial guess: Modify file coords and ncoord in the control-file accordingly, set riter to 0 and restart.

[prev] [prev-tail] [front] [up]