Documentation

docs/source/how-it-works.rst

@@ -55,6 +55,8 @@ The key function that any ``Engine`` subclass must have is ``calc_new(coords)``,
 .. note::
     The ``tests/test_customengine.py`` module shows how almost any energy function could be used to calculate the energy and gradient.  This was originally added by the `PySCF <https://pyscf.org/user/geomopt.html>`_ developers as an interface to their software.
 
+.. _internal_coordinate:
+
 Internal coordinate setup
 -------------------------
 
@@ -98,6 +100,8 @@ The diagonal elements are assigned following `Schlegel et al <https://link.sprin
 The initial diagonal elements for translation and rotation coordinates are set to 0.05.
 On the other hand, the default behavior in :ref:`transition state <transition>` optimizations is to calculate the full Hessian by finite difference of the gradient, because it is important to have the correct local shape of the PES for those calculations.
 
+.. _optimization_step:
+
 Optimization step
 -----------------
 
@@ -160,6 +164,8 @@ The iterations are repeated until the IC change from the initial coordinates mat
 
 The Cartesian step that matches the IC is obtained as :math:`\boldsymbol{\delta}^{(x)} = \mathbf{x}_n - \mathbf{x}_0`.
 
+.. _optimization_stepsize:
+
 Controlling the step size
 """""""""""""""""""""""""
 
@@ -220,6 +226,8 @@ The step size control algorithm can be summarized as:
 4. Iterate through Brent's method until :math:`f(r_{IC})` is converged to within 10% of the trust radius.
 5. Use the converged :math:`\boldsymbol{\delta}^{(x)}` to update the Cartesian coordinates for the next energy and gradient calculation.
 
+.. _convergence_criteria:
+
 Convergence criteria
 --------------------
 
@@ -254,13 +262,15 @@ The convergence criteria may be changed either as a set, by passing ``--converge
 Additionally, geomeTRIC supports Q-Chem-style and Molpro-style convergence criteria following the conventions of these packages. 
 Q-Chem requires the RMS gradient and *either* the RMS displacement or the energy change to fall below their criteria.
 Similarly, Molpro requires the maximum gradient and *either* the maximum displacement or the energy change to fall below their criteria.
-These alternate convergence conditions can be activated by passing ``--qccnv yes`` or ``-molpro yes`` on the command line.
+These alternate convergence conditions can be activated by passing ``--qccnv yes`` or ``--molpro yes`` on the command line.
 
 Trust radius adjustment
 -----------------------
 
 If the calculation is not converged yet, the step quality is calculated as:
 
+.. _step_quality:
+
 .. math::
     \begin{aligned}
     & Q = \frac{\Delta E_{\mathrm{actual}}}{\Delta E_{\mathrm{pred}}} \\

docs/source/install.rst

@@ -88,3 +88,20 @@ The Work Queue library in the `CCTools <https://github.com/cooperative-computing
 Installation of ``cctools`` can be done easily via conda as follows ::
 
     conda install -c conda-forge ndcctools
+
+.. _installbigchem:
+
+Installation of BigChem
+-----------------------
+
+BigChem can be used to perform the Hessian and NEB calculations in parallel.
+The github repository of `BigChem <https://github.com/mtzgroup/bigchem>`_ has detailed documentation explaining how its backend and broker work.
+
+BigChem can be installed using ``pip``::
+
+    pip install bigchem
+
+or it can be installed along with geomeTRIC::
+
+    pip install geometric[bigchem]
+

docs/source/irc.rst

@@ -0,0 +1,126 @@
+.. _irc:
+
+Intrinsic reaction coordinate
+=============================
+
+Basics
+------
+
+The Intrinsic Reaction Coordinate (IRC) method aims to trace a minimum energy pathway on a potential energy surface (PES) starting with an optimized transition state (TS) structure.
+This optimized TS structure sits at a first-order saddle point on the PES where the structure has only one imaginary vibrational frequency mode.
+To begin, geomeTRIC calculates the Hessian and vibrational frequencies to confirm this imaginary mode. The calculation of Hessian can be carried in parallel using :ref:`BigChem <installbigchem>` or :ref:`Work Queue <installcctools>`.
+
+Once the Hessian calculation is completed, the first step is taken in the positive direction of the corresponding eigenvector of the imaginary mode using mass-weighted internal coordinates.
+The Hessian is updated using the :ref:`BFGS method <bfgs-update>` and the succeeding steps are guided by the instantaneous acceleration vectors.
+The IRC method continues taking the steps until it meets the same convergence criteria as :ref:`geometry optimization <convergence_criteria>`.
+
+geomeTRIC will repeat the same procedure for the negative direction upon the convergence of the positive direction IRC and concatenate the two paths.
+
+Usage
+-----
+
+The IRC method can be used by passing ``--irc yes`` on the command line with ``geometric-optimize``. The input geometry should be an optimize TS structure.
+geomeTRIC adjusts the IRC step size based on the step quality, and the maximum step size will be set equal to the initial trust radius ``--trust [0.3]``.
+Once the both directions converge, geomeTRIC will write an output xyz file `input_irc.xyz` that shows the IRC pathway.
+
+Example
+-------
+
+See ```examples/1-simple-examples/hcn_hnc_irc``` directory.
+
+
+Theory
+------
+
+The IRC method uses a large portion of the same code as the optimization algorithm. The two main differences are in obtaining steps and adjusting the trust radius.
+
+Obtaining the IRC step
+""""""""""""""""""""""
+
+geomeTRIC implemented Gonzalez and Schlegel's `mass-weighted internal coordinates IRC method <https://doi.org/10.1021/j100377a021>`_.
+The mass-weighted Wilson B-matrix (:math:`\mathbf{B}`) has elements of :math:`dq_i / (dx_j \sqrt{m_j})`, where :math:`m_j` is the atomic mass, and the G-matrix is calculated as :math:`\mathbf{G} = \mathbf{B}\mathbf{B}^T` (See :ref:`internal coordinate setup <internal_coordinates>`).
+
+At the beginning of the first iteration, the mass-weighted step size (:math:`\mathbf{s}`) is calculated from the trust radius (:math:`R_{\mathrm{trust}}`) as follows:
+
+.. math::
+    \begin{aligned}
+    A = \sqrt{\sum_{i=1}^{N_{\mathrm{atoms}}} (\Delta x_i^2 + \Delta y_i^2 + \Delta z_i^2) \times \sqrt{m_i}} \\
+    \mathbf{s} = R_{\mathrm{trust}} \times A
+    \end{aligned}
+
+where :math:`\Delta x`, :math:`\Delta y`, and :math:`\Delta z` are the normalized Cartesian displacements along the imaginary mode on the saddle-point.
+The conversion factor :math:`A` is used in every iteration to convert the trust radius.
+
+Each IRC step starts with taking a half-step towards a pivot point following the instantaneous accelerations (:math:`-\mathbf{G} \cdot \mathbf{g}`).
+The pivot point (:math:`\mathbf{\mathrm{q}}^*`) is obtained by:
+
+.. math::
+    \begin{aligned}
+    \mathbf{\mathrm{q}}^* = \mathbf{\mathrm{q}} - 1/2 \mathbf{s} N \mathbf{G} \cdot \mathbf{g} \\
+    N = (\mathbf{g}^T \cdot \mathbf{G} \cdot \mathbf{g})^{1/2}
+    \end{aligned}
+
+where :math:`\mathbf{g}` is the gradients of internal coordinates :math:`\mathbf{\mathrm{q}}`.
+To reach the next point on the reaction path, another half-step needs to be taken from the pivot point along a vector that is 1) parallel to the acceleration vector of the next point and 2) has a scalar value equal to half of the mass-weighted step size (:math:`1/2\mathbf{s}`).
+
+First, a guessed point (:math:`\mathbf{\mathrm{q}}_1^\prime`) is obtained by taking the same step as the initial half-step starting from the pivot point.
+The guessed point is guided to the next guessed point (:math:`\mathbf{\mathrm{q}}_2^\prime`) until the vector pointing from the pivot point to the guessed point satisfies the two conditions.
+
+If we define the following two vectors:
+
+.. math::
+    \begin{aligned}
+    \mathbf{\mathrm{p}} = \mathbf{\mathrm{q}}_n^\prime - \mathbf{\mathrm{q}}^* \\
+    \Delta\mathbf{\mathrm{q}} = \mathbf{\mathrm{q}}_{n+1}^\prime - \mathbf{\mathrm{q}}_n^\prime
+    \end{aligned}
+    :label: vectors
+
+:math:`\Delta\mathbf{\mathrm{q}}` can be updated while keeping the scaler value of :math:`\mathbf{\mathrm{p}}` equal to :math:`1/2\mathbf{s}` until :math:`\mathbf{\mathrm{q}}_n^\prime` reaches the point where :math:`\mathbf{\mathrm{p}} + \Delta\mathbf{\mathrm{q}}` is parallel to the acceleration vectors at point :math:`\mathbf{\mathrm{q}}_{n+1}^\prime`.
+The vectors, Hessian, and gradients in mass-weighted internal coordinates can be expressed as:
+
+.. math::
+    \begin{aligned}
+    \Delta\mathbf{\mathrm{q}}_\mathbf{\mathrm{M}} = \mathbf{G}^{-1/2} \Delta\mathbf{\mathrm{q}}\\
+    \mathbf{\mathrm{p}}_\mathbf{\mathrm{M}} = \mathbf{G}^{-1/2} \mathbf{\mathrm{p}}\\
+    \mathbf{\mathrm{g}}_\mathbf{\mathrm{M}} = \mathbf{G}^{1/2} \mathbf{\mathrm{g}}^{\prime}\\
+    \mathbf{\mathrm{H}}_\mathbf{\mathrm{M}} = \mathbf{G}^{1/2} \mathbf{\mathrm{H}} \mathbf{G}^{1/2}\\
+    \end{aligned}
+    :label: mwic
+
+where :math:`\mathbf{\mathrm{g}}^{\prime}` represents the estimated gradients at the point :math:`\mathbf{\mathrm{q}}_n^\prime`, using a quadratic expansion.
+:math:`\mathbf{G}` is calculated at :math:`\mathbf{\mathrm{q}}_n^\prime` as well.
+
+The step size constraint can be expressed as:
+
+.. math::
+    (\mathbf{\mathrm{p}}_\mathbf{\mathrm{M}} + \Delta\mathbf{\mathrm{q}}_\mathbf{\mathrm{M}})^{T}(\mathbf{\mathrm{p}}_\mathbf{\mathrm{M}} + \Delta\mathbf{\mathrm{q}}_\mathbf{\mathrm{M}}) = (1/2 \mathbf{s})^2
+    :label: const1
+
+The other condition is satisfied at the convergence point (the next point), when the following equation holds true:
+
+.. math::
+    (\mathbf{\mathrm{g}}_\mathbf{\mathrm{M}} - \lambda \mathbf{\mathrm{p}}_\mathbf{\mathrm{M}}) + (\mathbf{\mathrm{H}}_\mathbf{\mathrm{M}} - \lambda \mathbf{I})\Delta\mathbf{\mathrm{q}}_\mathbf{\mathrm{M}} = 0
+    :label: const2
+
+where :math:`\lambda` is the Lagrangian multiplier and :math:`\mathbf{I}` is the identity matrix.
+
+Eq. :eq:`const2` can be rearranged as follows:
+
+.. math::
+    \Delta\mathbf{\mathrm{q}}_\mathbf{\mathrm{M}} = -(\mathbf{\mathrm{H}}_\mathbf{\mathrm{M}} - \lambda \mathbf{I})^{-1}(\mathbf{\mathrm{g}}_\mathbf{\mathrm{M}} - \lambda \mathbf{\mathrm{p}}_\mathbf{\mathrm{M}})
+    :label: delqm
+
+:math:`\lambda` is calculated iteratively after introducing Eq. :eq:`delqm` to :eq:`const1`.
+:math:`\Delta\mathbf{\mathrm{q}}_\mathbf{\mathrm{M}}` is then used to move :math:`\mathbf{\mathrm{q}}_n^\prime` to :math:`\mathbf{\mathrm{q}}_{n+1}^\prime` and new Eq. :eq:`vectors` and Eq. :eq:`mwic` are defined to calculate the next :math:`\Delta\mathbf{\mathrm{q}}_\mathbf{\mathrm{M}}`.
+This process repeats until the norm of :math:`\Delta\mathbf{\mathrm{q}}` falls below 1e-6. It then takes the rest of the half-step along :math:`\mathbf{\mathrm{p}} + \Delta\mathbf{\mathrm{q}}` from the pivot point, which completes an iteration.
+
+
+Trust radius adjustment
+"""""""""""""""""""""""
+
+The step quality (:math:`Q`) is calculated in the same way as the :ref:`energy minimization step quality <step_quality>`.
+The trust radius is adjusted as follows:
+
+* :math:`Q \geq 0.75` : "Good" step, trust radius is increased by a factor of :math:`\sqrt{2}`, but not greater than the maximum.
+* :math:`0.75 > Q \geq 0.50` : "Okay" step, trust radius is unchanged.
+* :math:`Q < 0.50` : Step is rejected, trust radius is decreased by setting it to :math:`0.5 \times \mathrm{min}(R_{\mathrm{trust}}, \mathrm{RMSD})`, but not lower than the minimum