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).

tsopt:
 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
calc:
 type: xtb
 pal: 4
 charge: 0
 mult: 1
geom:
 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.

opt:
 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.
calc:
 type: dimer     # Dimer calculator wrapping PySCF at HF/3-21G level of theory
 calc:
  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.
geom:
 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.

class pysisyphus.tsoptimizers.TSHessianOptimizer.TSHessianOptimizer(geometry, roots=None, root=0, hessian_ref=None, rx_modes=None, prim_coord=None, rx_coords=None, hessian_init='calc', hessian_update='bofill', hessian_recalc_reset=True, max_micro_cycles=50, trust_radius=0.3, trust_max=0.5, augment_bonds=False, min_line_search=False, max_line_search=False, assert_neg_eigval=False, **kwargs)

Bases: HessianOptimizer

Optimizer to find first-order saddle points.

__init__(geometry, roots=None, root=0, hessian_ref=None, rx_modes=None, prim_coord=None, rx_coords=None, hessian_init='calc', hessian_update='bofill', hessian_recalc_reset=True, max_micro_cycles=50, trust_radius=0.3, trust_max=0.5, augment_bonds=False, min_line_search=False, max_line_search=False, assert_neg_eigval=False, **kwargs)

Baseclass for transition state optimizers utilizing Hessian information.

Several arguments expect a typed primitive or an iterable of typed primitives. A typed primitive is specified as (PrimType, int, int, ...), e.g., for a bond between atoms 0 and 1: (BOND, 0, 1) or for a bend between the atom triple 0, 1, 2 as (BEND, 0, 1, 2).

Parameters:

• geometry (Geometry) -- Geometry to be optimized.

• roots (Optional[List[int]], default: None) -- Indices of modes to maximize along, e.g., to optimize saddle points of 2nd order. Overrides 'root'.

• root (int, default: 0) -- Index of imaginary mode to maximize along. Shortcut for 'roots' with only one root.

• hessian_ref (Optional[str], default: None) -- Filename pointing to a pysisyphus HDF5 Hessian.

• rx_modes (iterable of (iterable of (typed_prim, phase_factor))) -- Select initial root(s) by overlap with a modes constructed from the given typed primitives with respective phase factors.

• prim_coord (typed_prim) -- Select initial root/imaginary mode by overlap with this internal coordinate. Shortcut for 'rx_modes' with only one internal coordinate.

• rx_coords (iterable of (typed_prim)) -- Construct imaginary mode comprising the given typed prims by modifying a model Hessian.

• hessian_init (Literal['calc', 'unit', 'fischer', 'lindh', 'simple', 'swart', 'xtb', 'xtb1', 'xtbff'], default: 'calc') -- Type of initial model Hessian.

• hessian_update (Literal['none', None, False, 'bfgs', 'damped_bfgs', 'flowchart', 'bofill', 'ts_bfgs', 'ts_bfgs_org', 'ts_bfgs_rev'], default: 'bofill') -- Type of Hessian update. Defaults to BFGS for minimizations and Bofill for saddle point searches.

• hessian_recalc_reset (bool, default: True) -- Whether to skip Hessian recalculation after reset. Undocumented.

• max_micro_cycles (int, default: 50) -- Maximum number of RS iterations.

• trust_radius (float, default: 0.3) -- Initial trust radius in whatever unit the optimization is carried out.

• trust_max (float, default: 0.5) -- Maximum trust radius.

• augment_bonds (bool, default: False) -- Try to derive additional streching coordinates from the imaginary mode.

• min_line_search (bool, default: False) -- Carry out line search along the imaginary mode.

• max_line_search (bool, default: False) -- Carry out line search in the subspace that is minimized.

• assert_neg_eigval (bool, default: False) -- Check for the existences for at least one significant negative eigenvalue. If enabled and no negative eigenvalue is present the optimization will be aborted.

• **kwargs -- Keyword arguments passed to the HessianOptimizer/Optimizer baseclass.

static do_line_search(e0, e1, g0, g1, prev_step, maximize, logger=None)

prepare_opt(*args, **kwargs)

property root

property roots

step_and_grad_from_line_search(energy, gradient_trans, eigvecs, min_indices, max_indices)

update_ts_mode(eigvals, eigvecs)

valid_updates = ('bofill', 'ts_bfgs', 'ts_bfgs_org', 'ts_bfgs_rev')

12.5. TS-Optimizer using hessian information

12.5.1. Restricted-Step Partitioned-RFO

class pysisyphus.tsoptimizers.RSPRFOptimizer.RSPRFOptimizer(geometry, roots=None, root=0, hessian_ref=None, rx_modes=None, prim_coord=None, rx_coords=None, hessian_init='calc', hessian_update='bofill', hessian_recalc_reset=True, max_micro_cycles=50, trust_radius=0.3, trust_max=0.5, augment_bonds=False, min_line_search=False, max_line_search=False, assert_neg_eigval=False, **kwargs)

Bases: TSHessianOptimizer

optimize()

12.5.2. Restricted-Step Image-Function-RFO

class pysisyphus.tsoptimizers.RSIRFOptimizer.RSIRFOptimizer(geometry, roots=None, root=0, hessian_ref=None, rx_modes=None, prim_coord=None, rx_coords=None, hessian_init='calc', hessian_update='bofill', hessian_recalc_reset=True, max_micro_cycles=50, trust_radius=0.3, trust_max=0.5, augment_bonds=False, min_line_search=False, max_line_search=False, assert_neg_eigval=False, **kwargs)

Bases: TSHessianOptimizer

optimize()

12.5.3. Trust-Region Image-Function Method

class pysisyphus.tsoptimizers.TRIM.TRIM(geometry, roots=None, root=0, hessian_ref=None, rx_modes=None, prim_coord=None, rx_coords=None, hessian_init='calc', hessian_update='bofill', hessian_recalc_reset=True, max_micro_cycles=50, trust_radius=0.3, trust_max=0.5, augment_bonds=False, min_line_search=False, max_line_search=False, assert_neg_eigval=False, **kwargs)

Bases: TSHessianOptimizer

optimize()

12.6. TS-Optimization without hessian information

12.6.1. Dimer method

class pysisyphus.calculators.Dimer.Dimer(calculator, *args, N_raw=None, length=0.0189, rotation_max_cycles=15, rotation_method='fourier', rotation_thresh=0.0001, rotation_tol=1, rotation_max_element=0.001, rotation_interpolate=True, rotation_disable=False, rotation_disable_pos_curv=True, rotation_remove_trans=True, trans_force_f_perp=True, bonds=None, N_hessian=None, bias_rotation=False, bias_translation=False, bias_gaussian_dot=0.1, seed=None, write_orientations=True, forward_hessian=True, **kwargs)

property C

Shortcut for the curvature.

property N

add_gaussian(atoms, center, N, height=0.1, std=0.0529, max_cycles=50, dot_ref=None)

property can_bias_f0

property can_bias_f1

property coords0

property coords1

curvature(f1, f2, N)

Curvature of the mode represented by the dimer.

direct_rotation(optimizer, prev_step)

do_dimer_rotations(rotation_thresh=None)

property energy0

property f0

property f1

property f1_bias

property f2

Never calculated explicitly, but estimated from f0 and f1.

fourier_rotation(optimizer, prev_step)

get_N_raw_from_hessian(h5_fn, root=0)

get_forces(atoms, coords)

Meant to be extended.

get_gaussian_energies(coords, sum_=True)

get_gaussian_forces(coords, sum_=True)

get_hessian(atoms, coords)

Get Hessian matrix. Fall back to numerical Hessian, if not overriden.

Preferrably, this method should provide an analytical Hessian.

remove_translation(displacement)

property rot_force

rotate_coords1(rad, theta)

Rotate dimer and produce new coords1.

set_N_raw(coords)

property should_bias_f0

May lead to calculation of f0 and/or f1 if present!

property should_bias_f1

May lead to calculation of f0 and/or f1 if present!

update_orientation(coords)

class pysisyphus.calculators.Dimer.Gaussian(height, center, std, N)

energy(R, height=None)

forces(R, height=None)

exception pysisyphus.calculators.Dimer.RotationConverged