.. _tutorial_irc: 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 :ref:`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: :math:`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 :download:`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 :download:`here ` .. raw:: html :file: files/irc/irc_tutorial.html 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` .. include:: qa.inc