Intrinsic reaction coordinate

In this tutorial we will show how to use MLatom to find the reaction path that follows the intrinsic reaction coordinate (IRC). IRC starts from the transition state and there is specific tutorial about how to locate transition state with MLatom.

Theory behind

What is IRC

IRC (intrinsic reaction coordinate) is the minimal energy path (MEP) that connects reactant, product, and transition state (TS), in mass weighted coordinates. We can use an additional parameter $s$ as the arc length along the reaction path and define the IRC path as $x(s)$. It’s the solution to the equation:

\(dx(s)/ds = -g(s)/|g(s)|\)

which indicates that the direction of each point on IRC follows the direction of force (definition of MEP). It can be solved with standard numerical integration techniques.

It’s obvious from this equation that TS is not well described, because its energy gradient is zero. Thus, in IRC, the tangent of the curve at TS is the so-called transition vector, i.e., the eigenvector of the TS Hessian corresponding to the negative eigenvalue.

Why to use IRC

  • check whether the TS connects to the correct reactant and product.

  • check if intermediate exists

Example with MLatom

In the following jupyter notebook, we will show how to use MLatom to perform IRC calculation starting from the transition state of Diels-Alder reaction between ethene and 1,3-butadiene. The molecular file in .json format can be downloaded here. It stores results from frequency calculations including frequencies, force constants, normal modes, etc., which can be parsed and loaded to standard molecule object in MLatom. The full jupyter notebook can be downloaded here

irc_tutorial

Choice of programs

program can be used to specify different programs for IRC calculation and program_kwargs to pass settings to each program. Currently, the available options are None (the defualt one that uses our built-in implementation), pysisyphus, gaussian and geometric. In general, the recommended program is pysisyphus, which provides various algorithms and settings for customization. If you have dependency problems with it, we also have native implementations of LQA in MLatom to handle most cases.

Built-in implementation

Pros and Cons

  • Pros:

    • Compatible to MLatom and no additional installation

  • Cons

    • Only LQA (Local Quadratic Approximation) is supported

Useful keywords for IRC

  • n_steps (int): Defaults to 10.

  • step_size (float): Defaults to 0.1 Angstrom in mass-deweighted coordinates.

  • init_length (float): Defaults to 0.1 Angstrom in mass-deweighted coordinates.

  • hess_est (bool, str): Defaults to bofill . Available options are False (hessian calculation in each step), sr1, psb, bofill

pysisyphus

Pysisyphus is a Python package that is designed for exploration of potential energy surfaces in ground- and excited states. It can be easily installed with pip install pysisyphus.

Pros and Cons

  • Pros

    • rich algorithms and settings for customization

  • Cons

    • dependency problems

Useful keywords for IRC

Below we list some commonly used keywords. For full list of keywords available, please check: https://github.com/eljost/pysisyphus/blob/master/pysisyphus/irc/IRC.py

  • algorithm (str) : Available options are euler, eulerpc, gs, imk, lqa, modekill, rk4. Defaults to eulerpc

  • step_length (float): Defaults to 0.1 Bohr in mass weighted coordinates.

  • displ (str): Available options are energy, length, and energy_cubic. Defaults to energy

  • max_cycles (int): Defaults to 125

  • rms_grad_thresh (float): Defaults to 1e-3

  • energy_thresh (float): Defaults to 1e-6

Gaussian

Pros and Cons

  • Pros

    • The produced .log file can be recognized by many visualization softwares (MLatom can also load it with ml.molecule.load('myfile.log',format='gaussian'))

  • Cons (with interface in MLatom)

    • Heavy IO

    • Unable to reuse results of frequency calculation from MLatom

Keywords for IRC

The available keywords for IRC can be found here: https://gaussian.com/irc/. In our implementation, CalcFC is set by default. In this way, the hessian calculations can be handled by MLatom.

To use them in MLatom, they can be passed to irc_kwargs as program_kwargs={'irc_keywords': ['LQA']}. In the Gaussian input file, it will be written as irc(CalcFC,LQA).

geomeTRIC

geomeTRIC is not the recommended choice for IRC. It doesn’t provide any options to choose direction (at the time of implementation). The algorithm it uses is very different from the common one for IRC. To propagate IRC trajectories with geomeTRIC, it must be upgraded to version >=1.1.0

Any questions or suggestions?

If you have further questions, criticism, and suggestions, we would be happy to receive them in English or Chinese via email, Slack (preferred), or QQ (135258964). See our contact page for more information on how to subscribe to updates and follow us on social media.