neb and irc documentations are updated. Author/contributor informatio…

…n is removed from scripts.

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/irc.rst

@@ -14,7 +14,7 @@ Once the Hessian calculation is completed, the first step is taken in the positi
 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.
+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
 -----
@@ -26,7 +26,7 @@ Once the both directions converge, geomeTRIC will write an output xyz file `inpu
 Example
 -------
 
-See ```examples/1-simple-examples/hcn_hnc_irc``` directory.
+See ``examples/1-simple-examples/hcn_hnc_irc`` directory.
 
 
 Theory
@@ -38,14 +38,14 @@ 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>`).
+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}} \\
-    \mathbf{s} = R_{\mathrm{trust}} \times A
+    & 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.
@@ -55,64 +55,61 @@ Each IRC step starts with taking a half-step towards a pivot point following the
 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}
+    \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{\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}`).
+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{\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.
+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{\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
+    & \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{\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:`\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{\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}\\
+    & \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{\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.
+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{\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
+    (\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{\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
+    (\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{\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}})
+    \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{\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.
+: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

docs/source/neb.rst

@@ -5,9 +5,8 @@ Nudged elastic band
 
 Basics
 ------
-
-The nudged elastic band (NEB) is a chain-of-states method used to optimize a sequence of molecular structures (images) that connect two energy basins on a potential energy surface (PES).
-This optimized sequence of images approximates the minimum energy pathway (MEP) on the PES.
+The Nudged Elastic Band (NEB) method is a chain-of-states approach used to identify the minimum energy pathway (MEP) by optimizing a series of molecular structures (images) that connect two energy basins on a potential energy surface (PES).
+The resulting optimized sequence of images approximates the MEP, with the highest-energy image serving as a starting structure for :ref:`transition state optimization <transition>` to precisely locate the first-order saddle point.
 
 There are two main components of forces that guide each image closer to MEP. These forces, which each image experiences, can be expressed as:
 
@@ -33,7 +32,7 @@ Usage
 -----
 
 To run the NEB calculation with geomeTRIC, you need two input files:a QC input file and an xyz file containing the input chain coordinates.
-The input xyz file should have more frames than the specified number of images argument (``--images [11]``).
+The input xyz file must contain at least the numer of images specified by the ``--images [11]`` argument.
 These two input files can be provided by running ``geometric-neb qc.input chain.xyz`` on the command line.
 The chain coordinates from the input xyz file will override the molecular geometry from the QC input file.
 
@@ -43,7 +42,7 @@ If ``--optep yes`` is passed, the two endpoints of the input chain will be optim
 During optimization, geomeTRIC writes an xyz file of an image that climbs up towards the first-order saddle point(``qc.tsClimb.xyz``) in the working directory.
 Once the NEB calculation converges, the climbing image can be used as a high-quality initial guess for the :ref:`transition state optimization <transition>`.
 
-Since the singlepoint energy/gradient calculations of each images for :math:`\mathbf{F}_{\mathrm{PES}}` are independent of each other, they can be carried out in parallel using :ref:`Work Queue <installcctools>`, :ref:`BigChem <installbigchem>`, or `QCFractal <http://docs.qcarchive.molssi.org/projects/QCFractal/en/stable/>`_.
+Since the single point energy/gradient calculations for each images in :math:`\mathbf{F}_{\mathrm{PES}}` are independent, they can be carried out as separate parallel tasks using tools such as :ref:`Work Queue <installcctools>`, :ref:`BigChem <installbigchem>`, or `QCFractal <http://docs.qcarchive.molssi.org/projects/QCFractal/en/stable/>`_.
 
 .. note::
     geomeTRIC only supports Cartesian coordinate system for the NEB method.
@@ -52,13 +51,13 @@ Since the singlepoint energy/gradient calculations of each images for :math:`\ma
 Example
 -------
 
-See ```examples/1-simple-examples/hcn_hnc_neb``` directory.
+See ``examples/1-simple-examples/hcn_hnc_neb`` directory.
 
 
 Theory
 ------
 
-The NEB method follows the same procedures for :ref:`optimization steps <optimization_step>` and :ref:`step size control <optimization_stepsize>` as geometry optimizations.
+The NEB method follows the similar procedures for :ref:`optimization steps <optimization_step>` and :ref:`step size control <optimization_stepsize>` as geometry optimizations.
 The details of how each force component is applied will be explained here, along with how the step quality is calculated.
 
 Force components
@@ -70,28 +69,28 @@ The perpendicular and parallel forces are obtained as following:
 
 .. math::
     \begin{aligned}
-    & \mathbf{F}_{\mathrm{PES}}^{\perp}(\mathbf{r}_i) = \mathbf{F}_{\mathrm{PES}}(\mathbf{r}_i) - \mathbf{F}_{\mathrm{PES}}(\mathbf{r}_i) \cdot \mathbf{\hat{\tau}}_i\\
-    & \mathbf{F}_{\mathrm{spring}}^{\parallel}(\Delta \mathbf{r}_{i+1,i}, \Delta \mathbf{r}_{i-1, i}) = k[(\mathbf{r}_{i+1} - \mathbf{r}_i) - (\mathbf{r}_i - \mathbf{r}_{i-1})] \cdot \mathbf{\hat{\tau}}_i \mathbf{\hat{\tau}}_i
+    & \mathbf{F}_{\mathrm{PES}}^{\perp}(\mathbf{r}_i) = \mathbf{F}_{\mathrm{PES}}(\mathbf{r}_i) - \mathbf{F}_{\mathrm{PES}}(\mathbf{r}_i) \cdot \hat{\mathbf{\tau}}_i\\
+    & \mathbf{F}_{\mathrm{spring}}^{\parallel}(\Delta \mathbf{r}_{i+1,i}, \Delta \mathbf{r}_{i-1, i}) = k([(\mathbf{r}_{i+1} - \mathbf{r}_i) - (\mathbf{r}_i - \mathbf{r}_{i-1})] \cdot \hat{\mathbf{\tau}}_i) \hat{\mathbf{\tau}}_i
     \end{aligned}
 
-The tangent vector (:math:`\mathbf{\hat{\tau}}_i`) is defined as:
+The tangent vector (:math:`\hat{\mathbf{\tau}}_i`) is defined as:
 
 .. math::
-    \mathbf{\hat{\tau}}_i=
+    \hat{\mathbf{\tau}}_i=
     \begin{cases}
-        \mathbf{\hat{\tau}}_i^+ = \mathbf{r}_{i+1} - \mathbf{r}_i& \textrm{if}\qquad E_{i+1} > E_{i} > E_{i-1}\\
-        \mathbf{\hat{\tau}}_i^- = \mathbf{r}_i - \mathbf{r}_{i-1}& \textrm{if}\qquad E_{i+1} < E_{i} < E_{i-1}
+        \hat{\mathbf{\tau}}_i^+ = \mathbf{r}_{i+1} - \mathbf{r}_i& \textrm{if}\qquad E_{i+1} > E_{i} > E_{i-1}\\
+        \hat{\mathbf{\tau}}_i^- = \mathbf{r}_i - \mathbf{r}_{i-1}& \textrm{if}\qquad E_{i+1} < E_{i} < E_{i-1}
     \end{cases}
 
 where :math:`E_{i}` is the energy of :math:`i`-th image.
 
 For the images located at extrema, the following tangent is applied:
 
 .. math::
-    \mathbf{\hat{\tau}}_i=
+    \hat{\mathbf{\tau}}_i=
     \begin{cases}
-        \mathbf{\hat{\tau}}_i^+ \Delta E_i^{\mathrm{max}} + \mathbf{\hat{\tau}}_i^- \Delta E_i^{\mathrm{min}}  & \textrm{if}\qquad E_{i+1} > E_{i-1}\\
-        \mathbf{\hat{\tau}}_i^+ \Delta E_i^{\mathrm{min}} + \mathbf{\hat{\tau}}_i^- \Delta E_i^{\mathrm{max}}  & \textrm{if}\qquad E_{i+1} < E_{i-1}
+        \hat{\mathbf{\tau}}_i^+ \Delta E_i^{\mathrm{max}} + \mathbf{\hat{\tau}}_i^- \Delta E_i^{\mathrm{min}}  & \textrm{if}\qquad E_{i+1} > E_{i-1}\\
+        \hat{\mathbf{\tau}}_i^+ \Delta E_i^{\mathrm{min}} + \mathbf{\hat{\tau}}_i^- \Delta E_i^{\mathrm{max}}  & \textrm{if}\qquad E_{i+1} < E_{i-1}
     \end{cases}
 
 where
@@ -105,7 +104,7 @@ During the optimization, when the maximum RMS gradient of the chain falls below
 The climbing image receives a newly defined force, which is:
 
 .. math::
-    \mathbf{F}_i = -\nabla E(\mathbf{r}_{i_{\mathrm{max}}}) + 2 \nabla E(\mathbf{r}_{i_{\mathrm{max}}}) \cdot \mathbf{\hat{\tau}}_{i_{\mathrm{max}}}\mathbf{\hat{\tau}}_{i_{\mathrm{max}}}
+    \mathbf{F}_i = -\nabla E(\mathbf{r}_{i_{\mathrm{max}}}) + 2 (\nabla E(\mathbf{r}_{i_{\mathrm{max}}}) \cdot \hat{\mathbf{\tau}}_{i_{\mathrm{max}}})\hat{\mathbf{\tau}}_{i_{\mathrm{max}}}
 
 Once both the average and maximum gradient of :math:`i`-th image satisfy the convergence criteria, which are 0.025 and 0.05 eV/Ang, respectively by default, the image is locked.
 The locked images won't be moved until their gradients exceed the convergence criteria and they are unlocked.
@@ -119,7 +118,7 @@ The NEB method assesses step quality through changes in band energy and gradient
 .. math::
     Q_E =
     & \begin{cases}
-    & \frac{2\Delta E_{\mathrm{pred}} - \Delta E_{\mathrm{actual}}}{\Delta E_{\mathrm{pred}}} & \textrm{if }\qquad \Delta E_{\mathrm{actual}}, \Delta E_{\mathrm{pred}} > 0 \textrm{and} \Delta E_{\mathrm{actual}} > \Delta E_{\mathrm{pred}}\\
+    & \frac{2\Delta E_{\mathrm{pred}} - \Delta E_{\mathrm{actual}}}{\Delta E_{\mathrm{pred}}} & \textrm{if }\Delta E_{\mathrm{actual}} > \Delta E_{\mathrm{pred}} > 0\\
     & \frac{\Delta E_{\mathrm{actual}}}{\Delta E_{\mathrm{pred}}} & \textrm{else }
     & \end{cases} \\
 

geometric/engine.py

@@ -1,12 +1,9 @@
 """
 engine.py: Communicate with QM or MM software packages for energy/gradient info
 
-Copyright 2016-2020 Regents of the University of California and the Authors
+This code is part of geomeTRIC.
 
-Authors: Lee-Ping Wang, Chenchen Song
-
-Contributors: Yudong Qiu, Daniel G. A. Smith, Sebastian Lee, Qiming Sun,
-Alberto Gobbi, Josh Horton
+Copyright 2016-2024 Regents of the University of California and the Authors
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:

geometric/errors.py

@@ -2,13 +2,11 @@
 errors.py : errors module with definition of Customized Errors
 
 The classes and subclasses in this module defines a "tree" relation of exceptions,
-that can be used thoughout the code for a consistent error handling pattern.
+that can be used throughout the code for a consistent error handling pattern.
 
-Copyright 2016-2020 Regents of the University of California and the Authors
+This code is part of geomeTRIC.
 
-Authors: Lee-Ping Wang, Chenchen Song
-
-Contributors: Yudong Qiu
+Copyright 2016-2024 Regents of the University of California and the Authors
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:

geometric/info.py

@@ -2,13 +2,11 @@
 info.py : Contains logo, citation messages, and other general information
 
 The classes and subclasses in this module defines a "tree" relation of exceptions,
-that can be used thoughout the code for a consistent error handling pattern.
+that can be used throughout the code for a consistent error handling pattern.
 
-Copyright 2016-2020 Regents of the University of California and the Authors
+This code is part of geomeTRIC.
 
-Authors: Lee-Ping Wang, Chenchen Song
-
-Contributors:
+Copyright 2016-2024 Regents of the University of California and the Authors
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:

geometric/internal.py

@@ -1,11 +1,9 @@
 """
 internal.py: Internal coordinate systems
 
-Copyright 2016-2020 Regents of the University of California and the Authors
+This code is part of geomeTRIC.
 
-Authors: Lee-Ping Wang, Chenchen Song
-
-Contributors: 
+Copyright 2016-2024 Regents of the University of California and the Authors
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met: