13. Intrinsic Reaction Coordinate (IRC)

By default two paths are integrated in plus- and minus-direction of the imaginary mode (transition vector, TV) at the transition state (TS). If the IRC is started from a non-stationary point with non-vanishing gradient the downhill: True argument can be set to integrate only one path and skip the initial hessian calculation.

IRCs are integrated in mass-weighted coordinates and the default integrator is EulerPC.

The initial displacement from the TS is done by requiring a certain energy lowering \(\mathrm{d}E\) from moving along the TV and calculating the corresponding step length from a quadratic potential: \(\mathrm{d}E = \frac{1}{2} (k \cdot \mathrm{d}q^2)\), with \(k\) being the eigenvalue of the TV (imaginary mode) and \(\mathrm{d}q\) the required step length. The default required energy lowering is \(0.0005 ~ \mathrm{au}\). Alternatively the initial can be done by a prescribed length: displ: length and displ_length: [desired step] (default is \(0.1 \sqrt{\mathrm{amu}} \cdot \mathrm{bohr}\)).

The resulting endpoints of the IRC integration can be further optimized to stationary poins by adding the endopt: section (vide infra). By setting fragments: True in endopt separate fragments will be detected and optimized individually. This may be useful if the molecules are only weakly bound. By setting fragments: total the total system, as well as the separate fragments will be optimized.

By default IRC path data is dumped to dump_fn: irc_data.h5 every dump_every: 5 cycles. IRC paths can be plotted with pysisplot --irc.

13.1. YAML example(s)

Below you can find an example YAML-input including the most important options that the user may want to modify.

geom:
 fn: hfabstraction_ts_opt_xtb.xyz   # Input coordinates
calc:
 type: xtb                          # extended tight-binding calculator
 pal: 4
 charge: 0
 mult: 1
irc:
 type: eulerpc                      # Similar to EulerPC from Gaussian

 #rms_grad_thresh: 0.001            # Convergence threshold
 #displ: energy|length|energy_cubic # How to do the initial displacement
 #displ_energy: 0.001               # Energy lowering in au (Hartree)
 #displ_length: 0.1                 # Step length along the TV
 #forward: True
 #backward: True
 #downhill: False                   # Only integrate downhill, disables forward/backward
 #hessian_init: null                # Path to HDF5 Hessian file
 #displ_third_h5: null              # Path to HDF5 file containing third derivative data
endopt:
 #fragments: False|True|total       # Detect & optimize fragments separately. Default is
                                    # False. When set to 'total' the total system as well
                                    # as the fragments are optimized.
 do_hess: False                     # Frequency calculation at the end

Further examples for IRC calculations from .yaml input can be found here.

13.2. IRC base class

Base class for IRC integrators from which actual IRC integrators are derived.

class pysisyphus.irc.IRC.IRC(geometry, step_length=0.1, max_cycles=125, downhill=False, forward=True, backward=True, root=0, hessian_init=None, displ='energy', displ_energy=0.001, displ_length=0.1, displ_third_h5=None, rms_grad_thresh=0.001, hard_rms_grad_thresh=None, energy_thresh=1e-06, imag_below=0.0, force_inflection=True, check_bonds=True, out_dir='.', prefix='', dump_fn='irc_data.h5', dump_every=5)

__init__(geometry, step_length=0.1, max_cycles=125, downhill=False, forward=True, backward=True, root=0, hessian_init=None, displ='energy', displ_energy=0.001, displ_length=0.1, displ_third_h5=None, rms_grad_thresh=0.001, hard_rms_grad_thresh=None, energy_thresh=1e-06, imag_below=0.0, force_inflection=True, check_bonds=True, out_dir='.', prefix='', dump_fn='irc_data.h5', dump_every=5)

Base class for IRC calculations.

Parameters:

• geometry (Geometry) -- Transtion state geometry, or initial geometry for downhill run.

• step_length (float, optional) -- Step length in unweighted coordinates.

• max_cycles (int, optional) -- Positive integer, controlloing the maximum number of IRC steps taken in a direction (forward/backward/downhill).

• downhill (bool, default=False) -- Downhill run from a non-stationary point with non-vanishing gradient. Disables forward and backward runs.

• forward (bool, default=True) -- Integrate IRC in positive s direction.

• backward (bool, default=True) -- Integrate IRC in negative s direction.

• root (int, default=0) -- Use n-th root for initial displacement from TS.

• hessian_init (str, default=None) -- Path to Hessian HDF5 file, e.g., from a previous TS calculation.

• displ (str, one of ("energy", "length", "energy_cubic")) -- Controlls initial displacement from the TS. 'energy' assumes a quadratic model, from which a step length for a given energy lowering (see 'displ_energy') is determined. 'length' corresponds to a displacement along the transition vector. 'energy_cubic' considers 3rd derivatives of the energy along the transition vector.

• displ_energy (float, default=1e-3) -- Required energy lowering from the TS in au (Hartree). Used with 'displ: energy|energy_cubic'.

• displ_length (float, default=0.1) -- Step length along the transition vector. Used only with 'displ: length'.

• displ_third_h5 (str, optional) -- Path to HDF5 file containing 3rd derivative information. Used with 'displ: energy_cubic'.

• rms_grad_thresh (float, default=1e-3,) -- Convergence is signalled when to root mean square of the unweighted gradient is less than or equal to this value

• energy_thresh (float, default=1e-6,) -- Signal convergence when the energy difference between two points is equal to or less than 'energy_thresh'.

• imag_below (float, default=0.0) -- Require the wavenumber of the imaginary mode to be below the given threshold. If given, it should be a negative number.

• force_inflection (bool, optional) -- Don't indicate convergence before passing an inflection point.

• check_bonds (bool, optional, default=True) -- Report whether bonds are formed/broken along the IRC, w.r.t the TS.

• out_dir (str, optional) -- Dump everything into 'out_dir' directory instead of the CWD.

• prefix (str, optional) -- Short string that is prepended to all files that are created by this class, e.g., trajectories and HDF5 dumps.

• dump_fn (str, optional) -- Base name for the HDF5 files.

• dump_every (int, optional) -- Dump to HDF5 every n-th cycle.

property coords

dump_data(dump_fn=None, full=False)

dump_ends(path, prefix, coords=None, trj=False)

property energy

get_conv_fact(mw_grad, min_fact=2.0)

get_endpoint_and_ts_geoms()

get_full_irc_data()

get_irc_data()

get_path_for_fn(fn)

property gradient

initial_displacement()

Returns non-mass-weighted steps in +s and -s direction for initial displacement from the TS. Earlier version only returned one step, that was later multiplied by either 1 or -1, depending on the desired IRC direction (forward/backward). The current implementation directly returns two steps for forward and backward direction. Whereas for plus and minus steps for displ 'length' and displ 'energy'

step_plus = -step_minus

is valid, it is not valid for dipsl 'energy_cubic' anymore. The latter step is formed as

x(ds) = ds * v0 + ds**2 * v1

so

x(ds) != -x(ds)

as

ds * v0 + ds**2 * v1 != -ds * v0 - ds**2 * v1 .

So, all required step are formed directly and later used as appropriate.

See

https://aip.scitation.org/doi/pdf/10.1063/1.454172 https://pubs.acs.org/doi/10.1021/j100338a027 https://aip.scitation.org/doi/pdf/10.1063/1.459634

irc(direction)

log(msg)

property m_sqrt

mass_weigh_hessian(hessian)

property mw_coords

property mw_gradient

postprocess()

prepare(direction)

report_bonds(prefix, bonds)

report_conv_thresholds()

run()

set_data(prefix)

unweight_vec(vec)

valid_displs = ('energy', 'length', 'energy_cubic')

13.3. IRC Integrators

13.3.1. Damped-Velocity-Verlet integrator

class pysisyphus.irc.DampedVelocityVerlet.DampedVelocityVerlet(geometry, v0=0.04, dt0=0.5, error_tol=0.003, max_cycles=150, **kwargs)

Bases: IRC

damp_velocity(velocity)

estimate_error(new_mw_coords)

mw_grad_to_acc(mw_grad)

Takes care of the units for the mass-weighted gradient. Converts units of a mass-weighted gradient [Hartree/(Bohr*amu)] to units of acceleration [sqrt(amu)*Bohr/fs²]. The 1e30 comes from converting second² to femto second².

prepare(direction)

step()

13.3.2. Euler integrator

Not recommended as it only produces reasonable results with very small step sizes.

class pysisyphus.irc.Euler.Euler(geometry, step_length=0.01, **kwargs)

Bases: IRC

step()

13.3.3. Euler-Predictor-Corrector integrator

Recommended IRC integrator and default choice.

class pysisyphus.irc.EulerPC.EulerPC(*args, hessian_recalc=None, hessian_update='bofill', hessian_init='calc', max_pred_steps=500, loose_cycles=3, dump_dwi=False, scipy_method=None, corr_func='mbs', **kwargs)

Bases: IRC

corrector_step(init_mw_coords, step_length, dwi)

get_integration_length_func(init_mw_coords)

prepare(*args, **kwargs)

scipy_corrector_step(init_mw_coords, step_length, dwi)

Solve IRC equation dx/ds = -g/|g| on DWI PES in mass-weighted coordinates. Integration done until self.step_length in unweighted coordinates is achieved.

step()

13.3.4. Gonzales-Schlegel-2 integrator

class pysisyphus.irc.GonzalezSchlegel.GonzalezSchlegel(geometry, max_micro_cycles=20, micro_step_thresh=0.001, hessian_recalc=None, line_search=False, **kwargs)

Bases: IRC

micro_step(counter)

Constrained optimization on a hypersphere.

perp_component(vec, perp_to)

postprocess()

step()

13.3.5. Local-Quadratic-Approximation integrator

class pysisyphus.irc.LQA.LQA(geometry, N_euler=5000, **kwargs)

Bases: IRC

step()

13.3.6. Modified-Ishida-Morokuma-Komornicki integrator

Similar to the algorithm implemented by ORCA 4.

class pysisyphus.irc.IMKMod.IMKMod(geometry, corr_first=True, corr_first_size=0.5, corr_first_energy=True, corr_bisec_size=0.0025, corr_bisec_energy=True, **kwargs)

Bases: IRC

fit_grad_and_energies(energy0, grad0, energy1, direction)

fit_parabola(x, y)

step()

13.3.7. Runge-Kutta-4 integrator

Not recommended, as it is very slow.

class pysisyphus.irc.RK4.RK4(geometry, step_length=0.1, max_cycles=125, downhill=False, forward=True, backward=True, root=0, hessian_init=None, displ='energy', displ_energy=0.001, displ_length=0.1, displ_third_h5=None, rms_grad_thresh=0.001, hard_rms_grad_thresh=None, energy_thresh=1e-06, imag_below=0.0, force_inflection=True, check_bonds=True, out_dir='.', prefix='', dump_fn='irc_data.h5', dump_every=5)

Bases: IRC

get_k(mw_coords)

step()