Q-Chem 5.0 User’s Manual

9.2 Improved Algorithms for Transition-Structure Optimization

Transition-structure searches tend to be more difficult (meaning, more likely to be unsuccessful) as compared to minimum-energy (equilibrium) geometry optimizations. Odds of success can be enhanced via an initial guess structure that is determined in an automated way, rather than simply “guessed” by the user. Several such automated algorithms are available in Q-Chem, and are described in this section.

9.2.1 Freezing String Method

Perhaps the most significant difficulty in locating transition states is to obtain a good initial guess of the geometry to feed into a surface-walking algorithm. This difficulty becomes especially relevant for large systems, for which the dimensionality of the search space is large. Interpolation algorithms are promising for locating good guesses of the minimum-energy pathway connecting reactant and product states as well as approximate saddle-point geometries. For example, the nudged elastic band method [536, 537] and the string method [538] start from a certain initial reaction pathway connecting the reactant and the product state, and then optimize in discretized path space towards the minimum-energy pathway. The highest-energy point on the approximate minimum-energy pathway becomes a good initial guess for the saddle-point configuration that can subsequently be used with any local surface-walking algorithm.

Inevitably, the performance of any interpolation method heavily relies on the choice of the initial reaction pathway, and a poorly-chosen initial pathway can cause slow convergence, or possibly convergence to an incorrect pathway. The growing-string method [539] and freezing-string method [540, 541] offer solutions to this problem, in which two string fragments (one representing the reactant state and the other representing the product state) are “grown” (i.e., increasingly-finely defined) until the two fragments join. The freezing-string method offers a choice between Cartesian interpolation and linear synchronous transit (LST) interpolation. It also allows the user to choose between conjugate gradient and quasi-Newton optimization techniques.

Freezing-string calculations are requested by setting JOBTYPE = FSM in the $rem section. Additional job-control keywords are described below, along with examples. Consult Refs. Behn:2011 and Sharada:2012 for a guide to a typical use of this method.

FSM_NNODE

Specifies the number of nodes along the string


TYPE:

INTEGER


DEFAULT:

Undefined


OPTIONS:

$N$

number of nodes in FSM calculation


RECOMMENDATION:

$N=15$. Use 10 to 20 nodes for a typical calculation. Reaction paths that connect multiple elementary steps should be separated into individual elementary steps, and one FSM job run for each pair of intermediates. Use a higher number when the FSM is followed by an approximate-Hessian based transition state search (Section 9.2.2).


FSM_NGRAD

Specifies the number of perpendicular gradient steps used to optimize each node


TYPE:

INTEGER


DEFAULT:

Undefined


OPTIONS:

$N$

Number of perpendicular gradients per node


RECOMMENDATION:

Anything between 2 and 6 should work, where increasing the number is only needed for difficult reaction paths.


FSM_MODE

Specifies the method of interpolation


TYPE:

INTEGER


DEFAULT:

2


OPTIONS:

1

Cartesian

2

LST


RECOMMENDATION:

In most cases, LST is superior to Cartesian interpolation.


FSM_OPT_MODE

Specifies the method of optimization


TYPE:

INTEGER


DEFAULT:

Undefined


OPTIONS:

1

Conjugate gradients

2

Quasi-Newton method with BFGS Hessian update


RECOMMENDATION:

The quasi-Newton method is more efficient when the number of nodes is high.


An example input appears below. Note that the $molecule section includes geometries for two optimized intermediates, separated by ****. The order of the atoms is important, as Q-Chem assumes that the $n$th atom in the reactant moves toward the $n$th atom in the product. The FSM string is printed out in the file stringfile.txt, which contains Cartesian coordinates of the structures that connect reactant to product. Each node along the path is labeled in this file, and its energy is provided. The geometries and energies are also printed at the end of the Q-Chem output file, where they are labeled:

----------------------------------------
   STRING
----------------------------------------

Finally, if MOLDEN_FORMAT is set to TRUE, then geometries along the string are printed in a MolDen-readable format at the end of the Q-Chem output file. The highest-energy node can be taken from this file and used to run a transition structure search as described in section 9.1. If the string returns a pathway that is unreasonable, check whether the atoms in the two input geometries are in the correct order.

Example 9.211  Example of the freezing-string method.

$molecule
   0  1
   Si   1.028032  -0.131573  -0.779689
   H    0.923921  -1.301934   0.201724
   H    1.294874   0.900609   0.318888
   H   -1.713989   0.300876  -0.226231
   H   -1.532839   0.232021   0.485307
****
   Si   0.000228  -0.000484  -0.000023
   H    0.644754  -1.336958  -0.064865
   H    1.047648   1.052717   0.062991
   H   -0.837028   0.205648  -1.211126
   H   -0.855603   0.079077   1.213023
$end

$rem
   JOBTYPE        fsm
   FSM_NGRAD      3
   FSM_NNODE      12
   FSM_MODE       2
   FSM_OPT_MODE   2
   METHOD         b3lyp
   BASIS          6-31G
$end

9.2.2 Hessian-Free Transition-State Search

Once a guess structure to the transition state is obtained, standard eigenvector-following methods such as Baker’s partitioned rational-function optimization (P-RFO) algorithm [542] can be employed to refine the guess to the exact transition state. The reliability of P-RFO depends on the quality of the Hessian input, which enables the method to distinguish between the reaction coordinate (characterized by a negative eigenvalue) and the remaining degrees of freedom. In routine calculations therefore, an exact Hessian is determined via frequency calculation prior to the P-RFO search. Since the cost of evaluating an exact Hessian typically scales one power of system size higher than the energy or the gradient, this step becomes impractical for systems containing large number of atoms.

The exact Hessian calculation can be avoided by constructing an approximate Hessian based on the output of FSM.[535] The tangent direction at the transition state guess on the FSM string is a good approximation to the Hessian eigenvector corresponding to the reaction coordinate. The tangent is therefore used to calculate the correct eigenvalue and corresponding eigenvector by variationally minimizing the Rayleigh-Ritz ratio [543]. The reaction coordinate information is then incorporated into a guess matrix which, in turn, is obtained by transforming a diagonal matrix in delocalized internal coordinates [544, 545] to Cartesian coordinates. The resulting approximate Hessian, by design, has a single negative eigenvalue corresponding to the reaction coordinate. This matrix is then used in place of the exact Hessian as input to the P-RFO method.

An example of this one-shot, Hessian-free approach that combines the FSM and P-RFO methods in order to determine the exact transition state from reactant and product structures is shown below:

Example 9.212 

$molecule
   0  1
   Si   1.028032  -0.131573  -0.779689
   H    0.923921  -1.301934   0.201724
   H    1.294874   0.900609   0.318888
   H   -1.713989   0.300876  -0.226231
   H   -1.532839   0.232021   0.485307
****
   Si   0.000228  -0.000484  -0.000023
   H    0.644754  -1.336958  -0.064865
   H    1.047648   1.052717   0.062991
   H   -0.837028   0.205648  -1.211126
   H   -0.855603   0.079077   1.213023
$end

$rem
   JOBTYPE       fsm
   FSM_NGRAD     3
   FSM_NNODE     18
   FSM_MODE      2
   FSM_OPT_MODE  2
   METHOD        b3lyp
   BASIS         6-31g
   SYMMETRY      false
   SYM_IGNORE    true
$end

@@@

$rem
   JOBTYPE               ts
   SCF_GUESS             read
   GEOM_OPT_HESSIAN      read
   MAX_SCF_CYCLES        250
   GEOM_OPT_DMAX         50
   GEOM_OPT_MAX_CYCLES   100
   METHOD                b3lyp
   BASIS                 6-31g
   SYMMETRY              false
   SYM_IGNORE            true
$end

$molecule
   read
$end

9.2.3 Improved Dimer Method

Once a good approximation to the minimum energy pathway is obtained, e.g., with the help of an interpolation algorithm such as the growing string method, local surface walking algorithms can be used to determine the exact location of the saddle point. Baker’s P-RFO method [542], using either an approximate or an exact Hessian, has proven to be a very powerful for this purpose, but does require calculation of a full Hessian matrix.

The dimer method [546], on the other hand, is a mode-following algorithm that requires only the curvature along one direction in configuration space, rather than the full Hessian, which can be accomplished using only gradient evaluations. This method is thus especially attractive for large systems where a full Hessian calculation might be prohibitively expensive, or for saddle-point searches where the initial guess is such that the eigenvector of corresponding to the smallest Hessian eigenvalue does not correspond to the desired reaction coordinate. An improved version of the original dimer method [547, 548] has been implemented in Q-Chem, which significantly reduces the influence of numerical noise and thus significantly reduces the cost of the algorithm.