Merge pull request #201 from hjnpark/doc

Documentation

README.md

@@ -8,7 +8,7 @@
 This is a geometry optimization code for molecular structures.
 The code works by calling external software for the energy and 
 gradient through wrapper functions.  Q-Chem, TeraChem, Psi4, 
-Molpro, and Gaussian 09/16 are supported quantum chemistry
+Molpro, Gaussian 09/16, and CFOUR are supported quantum chemistry
 codes through the command line interface.  The PySCF and
 QCArchive packages also provide interfaces to geomeTRIC for
 optimization. MM optimizations using OpenMM and Gromacs are

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
 """""""""""""""""""""""""
 
@@ -187,10 +193,10 @@ This is designed to minimize the number of times the IC step is converted to Car
 .. image:: images/taxol-step-test.png
    :width: 600
 
-The above image shows the relationship between :math:`r_{IC}` (x-axis) and the RMSD (left y-axis; blue) for the structure shown (the taxol example at the 10th optimization cycle)
+The above image shows the relationship between :math:`r_{IC}` (x-axis) and the RMSD (left y-axis; blue) for the structure shown (the taxol example at the 10th optimization cycle).
 The plot is generated by scanning the value of :math:`r_{IC}`, optimizing :math:`\lambda` for each value, the calculating :math:`\boldsymbol{\delta}^{(x)}` and the RMSD.
 
-The optimization of :math:`\lambda` for a target step size (i.e. :math:`|\boldsymbol{\delta}^{(q)}| \rightarrow r_{IC}` follows `Hebden et al. <https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.716.2997>`_
+The optimization of :math:`\lambda` for a target step size (i.e. :math:`|\boldsymbol{\delta}^{(q)}| \rightarrow r_{IC}`) follows `Hebden et al. <https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.716.2997>`_
 The derivative of the step size with respect to :math:`\lambda_n` at iteration :math:`n` is given by:
 
 .. math::
@@ -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/index.rst

@@ -26,6 +26,8 @@ Main features of geomeTRIC include:
    how-it-works
    constraints
    transition
+   neb
+   irc
    meci
    engines
    help

docs/source/install.rst

@@ -41,7 +41,7 @@ Dependencies
 
 The required packages for geomeTRIC are as follows. Older versions of packages may work, no guarantees.
 
-* Python : Versions 2.7, 3.5, 3.6, 3.7 are supported. Basic functionalities tested on 3.6 to 3.10.
+* Python : Versions 2.7, 3.5+ are supported. Basic functionalities tested on 3.7 to 3.13.
 * NumPy: Version 1.15 or above
 * NetworkX : Version 2.2 or above
 
@@ -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,123 @@
+.. _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 repeats the same procedure for the negative direction after the positive direction IRC converges. Once both directions have converged, one of the paths is reversed to concatenated with the other, using the TS structure as the connecting point.
+
+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_coordinate>`).
+
+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}} \\
+    & \mathrm{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::
+    \mathbf{q}^* = \mathbf{q} - \frac{\mathrm{s}}{2} \cdot \frac{\mathbf{G} \cdot \mathbf{g}}{(\mathbf{g}^T \cdot \mathbf{G} \cdot \mathbf{g})^{2}}
+
+where :math:`\mathbf{g}` is the gradients of internal coordinates :math:`\mathbf{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:`\mathrm{s}/2`).
+
+First, a guessed point (:math:`\mathbf{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{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{p} = \mathbf{q}_n^\prime - \mathbf{q}^* \\
+    & \Delta\mathbf{q} = \mathbf{q}_{n+1}^\prime - \mathbf{q}_n^\prime
+    \end{aligned}
+    :label: vectors
+
+:math:`\Delta\mathbf{q}` can be updated while keeping the scaler value of :math:`\mathbf{p}` equal to :math:`\mathrm{s}/2` until :math:`\mathbf{q}_n^\prime` reaches the point where :math:`\mathbf{p} + \Delta\mathbf{q}` is parallel to the acceleration vectors at point :math:`\mathbf{q}_{n+1}^\prime`.
+The vectors, Hessian, and gradients in mass-weighted internal coordinates can be expressed as
+
+.. math::
+    \begin{aligned}
+    & \Delta\mathbf{q}_\mathrm{M} = \mathbf{G}^{-1/2} \Delta\mathbf{q}\\
+    & \mathbf{p}_\mathrm{M} = \mathbf{G}^{-1/2} \mathbf{p}\\
+    & \mathbf{g}_\mathrm{M} = \mathbf{G}^{1/2} \mathbf{g}^{\prime}\\
+    & \mathbf{H}_\mathrm{M} = \mathbf{G}^{1/2} \mathbf{H} \mathbf{G}^{1/2}\\
+    \end{aligned}
+    :label: mwic
+
+where :math:`\mathbf{g}^{\prime}` represents the estimated gradients at the point :math:`\mathbf{q}_n^\prime`, using a quadratic expansion.
+:math:`\mathbf{G}` is calculated at :math:`\mathbf{q}_n^\prime` as well.
+
+The step size constraint can be expressed as:
+
+.. math::
+    (\mathbf{p}_\mathrm{M} + \Delta\mathbf{q}_\mathrm{M})^{T}(\mathbf{p}_\mathrm{M} + \Delta\mathbf{q}_\mathrm{M}) = (\frac{\mathrm{s}}{2})^2
+    :label: const1
+
+The other condition is satisfied at the convergence point (the next point), when the following equation holds true:
+
+.. math::
+    (\mathbf{g}_\mathrm{M} - \lambda \mathbf{p}_\mathrm{M}) + (\mathbf{H}_\mathrm{M} - \lambda \mathbf{I})\Delta\mathbf{q}_\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{q}_\mathrm{M} = -(\mathbf{H}_\mathrm{M} - \lambda \mathbf{I})^{-1}(\mathbf{g}_\mathrm{M} - \lambda \mathbf{p}_\mathrm{M})
+    :label: delqm
+
+:math:`\lambda` is calculated iteratively after introducing Eq. :eq:`delqm` to :eq:`const1`.
+:math:`\Delta\mathbf{q}_\mathrm{M}` is then used to move :math:`\mathbf{q}_n^\prime` to :math:`\mathbf{q}_{n+1}^\prime` and new Eq. :eq:`vectors` and Eq. :eq:`mwic` are defined to calculate the next :math:`\Delta\mathbf{q}_\mathrm{M}`.
+This process repeats until the norm of :math:`\Delta\mathbf{q}` falls below 1e-6. It then takes the rest of the half-step along :math:`\mathbf{p} + \Delta\mathbf{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