12. Transition State Optimization

To cite Frank Jensen: "Locating transition states (TSs) is black magic, especially in internal coordinates" and I can say this is true. The most promising TS optimizers in pysisyphus employ second derivative information (hessian) but locating TS is also possible using only first derivatives by means of the dimer method (DM).

TS optimizations should preferably be carried out in internal coordinates (type: redund|dlc). Before starting the actual TS search, one should always check the internal coordinates that pysisyphus sets up, by executing pysistrj [xyz] --internals, with [xyz] corresponding to the geometry you plan to use for the optimization. Check carefully if all supposed reaction coordinates are present. If they are missing define them manually with the add_prims key in the YAML input.

Calculation of the exact hessian can be avoided by using the DM by using rx_coords: [list of primitives] in combination with hessian_init: fischer|lindh|swart|simple. By using both options, a diagonal model hessian is calcualted and modified for use in a TS optimization.

The DM can be used with a powerful preconditioned LBFGS optimizer (type: lbfgs).

12.1. Hessian Based TS Optimization

12.1.1. YAML example

Below you can find an example YAML-input including the most important options that the user may want to modify for the Hessian-based optimizers (RS-I-RFO, RS-P-RFO, TRIM).

 type: rsirfo|rsprfo|trim         # Optimization algorithm

 do_hess: True                    # Calculate the hessian at the final geometry
                                  # after the optimization.

 hessian_recalc: 1                # Recalculate the exact hessian every n-th cylce
                                  # Very costly but a small number like 5 may be a good
                                  # idea.

 prim_coord: [BOND, 3, 18]        # Select the mode to follow uphill by overlap with
                                  # this primitive internal coordinate. Expects exactly
                                  # one typed primitive.

rx_modes: [[[[BOND, 3, 18], 1]]]  # Select initial mode based on overlaps with a
                                  # constructed mode(s). Can be seen as generalization of
                                  # prim_coord. Expects a a list of (lists of pairs, with each
                                  # pair comprising a typed primitive and a phase factor).
                                  # See examples/tsopt/{05..,06..} for examples.

 #rx_coords: [[BOND, 3, 18],]     # Obtain initial mode by modifying a model Hessian.  Has
                                  # to be used with 'hess_init: swart|fischer|lindh|xtb'
                                  # to avoid calculation of exact hessian.
                                  # Also requires 'redund' or 'dlc' to work.

 #root: 0                         # Follow the n-th imaginary mode uphill, minimize
                                  # along the others.

 #hessian_init: calc              # Initial hessian is calculated exactly.

 #hessian_update: bofill          # Bofill hessian-update. Should not be modified for
                                  # TS optimizations.

 #hessian_ref: [path]             # Expects a path to a precalculated hessian at any
                                  # level of theory. The emode to maximize along is selected
                                  # by highest overlap with imaginary mode from the reference
                                  # hessian.

 #max_micro_cycles: 50            # No. of micro cycles for the RS-variants. Does not apply
                                  # to TRIM.

 #trust_radius: 0.3               # Initial trust radius.
 #trust_max: 1.0                  # Max. trust radius
 #trust_min: 0.1                  # Min. trust radius
 type: xtb
 pal: 4
 charge: 0
 mult: 1
 type: redund
 fn: shaked.geom_000.xyz
 add_prims: [[24, 20], ]          # If using internal coordinates ALWAYS check the coordinates
                                  # that pysisyphus generates (pysistrj [xyz] --internals). If
                                  # some important (reaction) coordinates appears to be missing
                                  # define them manually. In this example
                                  # we add a bond. If additional internal coordinates can
                                  # derived from the added primitives pysisyphus will do
                                  # it.

Second-derivative-free TS optimization can be done using the Dimer method. Attached you can find a sample input.

Further examples for TS optimizations from .yaml input can be found here.

12.2. Dimer Method

pysisyphus implements the dimer method (DM) as calculator, wrapping another calculator for the actual energy and gradient calculations. All DM configurations have to be done in the calc: section. The best results are obtained by optimizing the DM with LBFGS in connection with a preconditioner (PLBFGS), estimated from the Lindh model Hessian. in constrast to the TS optimizers mentioned above, PBLFGS is configured in the opt section of the YAML input. DM related information is logged to dimer.log. The current DM orientation is saved in files with extension .N. Animated .trj files of the DM orientation are saved to .N.trj.

12.2.1. YAML example

Below you can find an example for the DM, applied to the isomerization of HCN. Default values are commented out. See examples/tsopt/02_hcn_tsopt_dimer for the full example.

 type: plbfgs    # Preconditioned LBFGS
 thresh: baker   # Baker threshold, comparable to the Gaussian defaults
 do_hess: True   # Calculate Hessian at the end. May not be applicable to systems
                 # where Hessian calculation is infeasible.
 type: dimer     # Dimer calculator wrapping PySCF at HF/3-21G level of theory
  type: pyscf
  basis: 321g
  pal: 2
  charge: 0
  mult: 1
 #N_raw: [file]                # Path to a file containing an initial dimer orientation
 #length: 0.0189               # Separation of dimer images from common midpoint
 #rotation_max_cycles: 15      # Maximum number of rotations
 #rotation_method: fourier
 #rotation_thresh: 1e-4        # rms(rotationa force) threshold
 #rotation_tol: 1              # Rotation tolerance in degree. If a proposed trial
                               # rotation angle is below the tolerance the rotation
                               # will be skipped.
 #rotation_max_element: 0.001  # Max. step element for "rotation_method: direct"
 #rotation_interpolate: True   # Interpolate force on image after rotation
 #rotation_remove_trans: True  # Remove overall translation from N-vector
 #seed: 20182503                # Seed for the RNG for reproducability reasons.
 type: cart
 fn: 01_hcn.xyz

12.3. General advice for TS optimizations

  • Use as many exact hessians as your computational budget allows (hessian_recalc: [n])

  • Use tight convergence settings for your SCF. In contrast to codes like ORCA pysisyphus does not enforces this automatically

  • Grid-based methods (DFT, some RI-types) may need finer grids to reduce numerical noise

  • Try to preoptimize the TS using a cheap method (xtb) that allows hessian_recalc: 1

  • Decrease current & allowed trust radius (trust_radius: 0.1 and trust_max: 0.3|0.5)

  • Check for orrect coordinate setup

12.4. TSHessianOptimizer base class

Base class for TS optimizers that empoly a hessian.

12.5. TS-Optimizer using hessian information

12.5.1. Restricted-Step Partitioned-RFO

12.5.2. Restricted-Step Image-Function-RFO

12.5.3. Trust-Region Image-Function Method

12.6. TS-Optimization without hessian information

12.6.1. Dimer method