diff --git a/docs/source/how-it-works.rst b/docs/source/how-it-works.rst index 186116b9..56fb0807 100644 --- a/docs/source/how-it-works.rst +++ b/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 `_ 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 ` 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}}} \\ diff --git a/docs/source/install.rst b/docs/source/install.rst index fd901589..069f70b1 100644 --- a/docs/source/install.rst +++ b/docs/source/install.rst @@ -88,3 +88,20 @@ The Work Queue library in the `CCTools `_ 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] + diff --git a/docs/source/irc.rst b/docs/source/irc.rst new file mode 100644 index 00000000..ad994772 --- /dev/null +++ b/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 ` or :ref:`Work Queue `. + +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 ` 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 `. + +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 `_. +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 `). + +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 `. +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 \ No newline at end of file diff --git a/docs/source/neb.rst b/docs/source/neb.rst new file mode 100644 index 00000000..31951222 --- /dev/null +++ b/docs/source/neb.rst @@ -0,0 +1,150 @@ +.. _neb: + +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. + +There are two main components of forces that guide each image closer to MEP. These forces, which each image experiences, can be expressed as: + +.. math:: + \mathbf{F}_i = \mathbf{F}_{\mathrm{PES}}^{\perp}(\mathbf{r}_i) + \mathbf{F}_{\mathrm{spring}}^{\parallel}(\Delta \mathbf{r}_{i+1,i}, \Delta \mathbf{r}_{i-1, i}) + :label: neb_force + +where :math:`\mathbf{r}_i` is the geometry of the :math:`i`-th image and :math:`\perp, \parallel` indicate that each force contribution is decomposed into perpendicular and parallel components to the path, and only the indicated component of each force contribution is included for NEB. +If all the components of both forces (:math:`\mathbf{F}_{\mathrm{PES}}` + :math:`\mathbf{F}_{\mathrm{spring}}`) are applied to the images, the chain is called a plain band. A hybrid band omits the parallel force from the PES (:math:`\mathbf{F}_{\mathrm{PES}}^{\perp}` + :math:`\mathbf{F}_{\mathrm{spring}}`). +geomeTRIC is capable of optimizing all three types of bands by minimizing the force components iteratively. + +The root-mean-squared (RMS) Cartesian gradients of each images (:math:`g_\mathrm{RMS}`) are calculated as: + + +.. math:: + g_\mathrm{RMS} = \sqrt{\frac{1}{N_{\mathrm{atoms}}} \displaystyle\sum_{i=1}^{N_{\mathrm{atoms}}} ({g_{x_i}}^2 + {g_{y_i}}^2 + {g_{z_i}}^2)} + :label: grms + +where :math:`g_{x_i}`, :math:`g_{y_i}`, and :math:`g_{z_i}` represent the gradients along the x, y, and z Cartesian axes, respectively. +The NEB calculation will continue iterating until the average and maximum :math:`g_\mathrm{RMS}` of all the images fall below convergence criteria. + +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]``). +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. + +By default, geomeTRIC will align all the images with the first image. This could be turned off by passing ``--align no``. +If ``--optep yes`` is passed, the two endpoints of the input chain will be optimized before alignment. + +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 `. + +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 `, :ref:`BigChem `, or `QCFractal `_. + +.. note:: + geomeTRIC only supports Cartesian coordinate system for the NEB method. + If ``--coordsys [tric]`` is specified and ``--optep yes`` is passed, the coordinate system will be used to optimize the two endpoints of the input chain. + +Example +------- + +See ```examples/1-simple-examples/hcn_hnc_neb``` directory. + + +Theory +------ + +The NEB method follows the same procedures for :ref:`optimization steps ` and :ref:`step size control ` 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 +"""""""""""""""" + +The gradients are calculated based on the force components described in Equation (:eq:`neb_force`). +geomeTRIC implemented Henkelman and Jónsson's `improved tangent `_ and `climibing image `_ method. +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 + \end{aligned} + +The tangent vector (:math:`\mathbf{\hat{\tau}}_i`) is defined as: + +.. math:: + \mathbf{\hat{\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} + \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= + \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} + \end{cases} + +where + +.. math:: + \Delta E_i^{\mathrm{max}} = max(|E_{i+1} - E_i|, |E_{i-1} - E_i|) \\ + \Delta E_i^{\mathrm{min}} = min(|E_{i+1} - E_i|, |E_{i-1} - E_i|) + +The tangent vector is normalized and applied to project the force components accordingly. +During the optimization, when the maximum RMS gradient of the chain falls below a threshold (default set at 0.5 ev/Ang using ``--climb [0.5]``), the highest energy image (:math:`i_{\mathrm{max}}`) is switched to climbing mode. +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}}} + +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. +The NEB calculation will converge when the average and maximum of RMS-gradient of all the images fall below the criteria. + +Trust radius adjustment +""""""""""""""""""""""" + +The NEB method assesses step quality through changes in band energy and gradients. The step quality based on energy (:math:`Q_E`) is expressed as: + +.. 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{\Delta E_{\mathrm{actual}}}{\Delta E_{\mathrm{pred}}} & \textrm{else } + & \end{cases} \\ + +where :math:`\Delta E_{\mathrm{actual}}` represents the energy difference between the previously iterated and current chain. +The :math:`\Delta E_{\mathrm{pred}}` is calculated using the same equation as the :ref:`predicted energy ` of the geometry optimization step. + +The step quality based on gradients (:math:`Q_g`) is calculated as: + +.. math:: + Q_g = 2 - \frac{g_{\mathrm{curr}}}{max(g_{\mathrm{pred}}, \frac{g_{\mathrm{prev}}}{2}, \frac{g_{\mathrm{conv}}}{2})} + +where :math:`g_{\mathrm{curr}}` and :math:`g_{\mathrm{prev}}` are the average RMS Cartesian gradients of current and previously iterated chain, respectively. +:math:`g_{\mathrm{conv}}` is the average gradient convergence criterion (``--avgg [0.025]``). +Predicted Cartesian gradients of each image (:math:`g_{\mathrm{cart}}`) is calculated as + +.. math:: + g_{\mathrm{cart}} = \mathbf{H} \cdot \boldsymbol \delta + \mathbf{g} + +where :math:`\boldsymbol \delta, \mathbf{g}, \mathbf{H}` are the displacement, gradient, and Hessian in Cartesian coordinate. +:math:`g_{\mathrm{pred}}` is the average RMS of :math:`g_{\mathrm{cart}}` calculated using Eq. :eq:`grms`. + +The larger value is chosen as the overall step quality (:math:`Q`) between the two values. +The overall step quality is then used to adjust the trust radius (:math:`R_{\mathrm{trust}}`) as following: + +* :math:`Q \geq 0.50` : "Good" step, trust radius is increased by a factor of :math:`\sqrt{2}`, but not greater than the maximum. +* :math:`0.50 > Q \geq 0.00` : "Okay" step, trust radius is unchanged. +* :math:`0.00 > Q \geq -1.0` : "Poor" step, 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. +* :math:`Q < -1.0` : Step is rejected in addition to decreasing the trust radius as above. diff --git a/docs/source/options.rst b/docs/source/options.rst index d2ed7758..d88f39fe 100644 --- a/docs/source/options.rst +++ b/docs/source/options.rst @@ -88,6 +88,14 @@ Also, a numerical Hessian will be computed at the start of the optimization (not .... +``--irc [yes/no]`` + +Provide ``yes`` to find an intrinsic reaction coordinate pathway starting from a first-order saddle point (optimized transition state structure). +It will compute a numerical Hessian at the beginning of the calculation. + +.... + + ``--rigid [yes/no]`` Provide ``yes`` to keep molecules / fragments rigid during the optimization. @@ -136,7 +144,7 @@ if the coordinates at run time matches the existing coordinate file, the Hessian Currently, Hessian matrices are computed by geomeTRIC by numerical central difference of the gradient, requiring 1+6*N(atoms) total calculations. The finite difference step size is ``1.0e-3`` a.u. which is the default in many QC programs. -Independent gradient jobs can be computed either locally in serial or in parallel using the Work Queue distributed computing library (see ``--port`` below). +Independent gradient jobs can be computed either locally in serial or in parallel using the Work Queue distributed computing library (see ``--wqport`` below). Individual gradient calculations are stored in folders such as ``[prefix].tmp/hessian/displace/001[m/p]`` which stands for "coordinate 1 minus/plus displacement". If the job is interrupted and restarted, existing completed gradient calculations will be read instead of recomputed. @@ -148,7 +156,7 @@ Interfaces for using native Hessian calculation routines will be added in the fu Possible values to pass to this argument are: - ``never`` (**default for minimization and MECI**) : Do not calculate the Hessian or read Hessian data. -- ``first`` (**default for transition state**) : Calculate the Hessian for the initial structure. +- ``first`` (**default for transition state and IRC**) : Calculate the Hessian for the initial structure. - ``last`` : Calculate the Hessian for the final (optimized) structure. - ``first+last`` : Calculate the Hessian for both the initial and final structure. - ``stop`` : Calculate the Hessian for the initial structure, and then stop (do not optimize). @@ -157,7 +165,14 @@ Possible values to pass to this argument are: .... -``--port [9876]`` +``--bigchem [yes/no]`` + +Provide ``yes`` to carry out the Hessian calculation in parallel using `BigChem `_. +Ensure that BigChem's backend and broker are running properly. + +.... + +``--wqport [9876]`` Provide a port number for the Work Queue distributed computing server. This is only used for distributing gradient calculations in numerical Hessian calculations. @@ -320,6 +335,113 @@ This is enabled by default in energy minimizations, and disabled in transition s If a number is provided, the internal coordinate system will be rebuilt at the specified interval as if the current structure were the input structure. This is disabled by default because it tends to lower performance, but may be useful for debugging. +.. _neb_options: + +NEB options +----------- + +These options control the NEB method. + +.... + +``--maxg [0.05]`` + +``--avgg [0.025]`` + +The convergence occurs when the maximum and average RMS-gradient of all images fall below these thresholds. + +.... + +``--guessk [0.05]`` + +This value will be used to build the guessed Hessian for the input chain. + +.... + +``--neb_history [1]`` + +Number of chain histories to save in memory. Note that each chain is memory intensive. + +.... + +``--neb_maxcyc [100]`` + +This sets the maximum number of NEB iterations. + +.... + +``--climb [0.5]`` + +``--ncimg [1]`` + +When the maximum RMS-gradient of the chain falls below ``--climb``, ``--ncimg`` number of images will be in climbing mode. +The climbing images will be selected in descending order from the highest energy image. + +.... + +``--plain [0]`` + +This option determines force components that will be projected. The default value is ``0`` which is NEB. +The other two available bands are hybrid (``1``) and plain (``2``) band. +The details of the force components of each band can be found in the :ref:` NEB page `. + +.... + +``--optep [yes/no]`` + +Provide ``yes`` to optimize the two end points of the input chain prior to the NEB calculation. It will not optimize the end points by default. + +.... + +``--align [yes/no]`` + +Provide ``yes`` to align all the images with the first image. It will align them by default. + +.... + +``--trust [0.1]`` + +``--tmax [0.3]`` + +These options control the starting and maximum values of the trust radius. +The trust radius is the maximum allowed Cartesian displacement of an optimization step measured in Angstrom. + +Depending on the quality of individual optimization steps, the trust radius can be increased from its current value up to the ``--tmax`` value, or it can be decreased down to a minimum value. + +The minimum trust radius cannot be user-set; its value is ``0.0`` for transition state and MECI jobs, and the smaller of the ``drms`` convergence criteria and ``1.2e-3`` for energy minimizations. +The purpose of the minimum trust radius is to prevent unpredictable behavior when the trust radius becomes extremely small (e.g. if the step is so small that the energy change is smaller than the SCF convergence criteria). + +.... + +``--skip [yes/no]`` + +Setting it to ``yes`` will skip Hessian updates that will introduce negative eigenvalues. + +.... + +``--epsilon [1e-5]`` + +When the eigenvalues of the Hessian fall below ``--epsilon``, it will rebuild the Hessian using the ``--guessk`` value. + +.... + +``--bigchem [yes/no]`` + +Provide ``yes`` to carry out the Hessian calculation in parallel using `BigChem `_. +Ensure that BigChem's backend and broker are running properly. + +.... + +``--wqport [9876]`` + +Provide a port number for the Work Queue distributed computing server. +This is only used for distributing gradient calculations in numerical Hessian calculations. +This number can range up to 65535, and a number in the high 4-digit range is acceptable. +Do not use privileged port numbers (less than 1024). + +The port number should not be used by other servers running on your machine, and should match +the port number provided to Work Queue worker processes whose job is to execute the gradient calculations. + Structure options ----------------- diff --git a/docs/source/transition.rst b/docs/source/transition.rst index 92d9e2eb..fe7a29f2 100644 --- a/docs/source/transition.rst +++ b/docs/source/transition.rst @@ -222,7 +222,7 @@ where i.e. :math:`\boldsymbol\delta` is the current optimization step, :math:`\boldsymbol\xi` is the difference between the actual gradient change and the *predicted* gradient change using the previous structure's Hessian, and :math:`\phi` measures the alignment between the two vectors, being equal to 1 when they are orthogonal and 0 when parallel. The PSB update is mixed in when the two vectors are almost orthogonal, as the MS update approaches zero and becomes unstable. -This updating method is used in TS optmization because it does not preserve positive-definiteness of the Hessian matrix, as opposed to BFGS (which is used in energy minimization). +This updating method is used in TS optimization because it does not preserve positive-definiteness of the Hessian matrix, as opposed to BFGS (which is used in energy minimization). All of the Hessian updates are carried out in internal coordinates. Step size control @@ -261,7 +261,7 @@ If you want to use an analytic Hessian from running the QC software separately, The Hessian file must be stored as a square matrix in Numpy-readable text format (not binary) with dimension :math:`3 \times N_{\mathrm{atoms}}`. The gradient calculations may be parallelized by distributing the jobs to remote "worker" nodes using the `Work Queue distributed computing library `_; this can greatly reduce the wall time relative to performing the gradient calculations serially. -To enable this behavior, first ensure that the Work Queue library and Python module are installed, then pass ``--port ####`` on the command line where ``####`` is a custom port number (I usually use a large four-digit number, such as 7953, that is not commonly used by other services). +To enable this behavior, first ensure that the Work Queue library and Python module are installed, then pass ``--wqport ####`` on the command line where ``####`` is a custom port number (I usually use a large four-digit number, such as 7953, that is not commonly used by other services). Then run the ``work_queue_worker`` program on the worker node, providing the host name that is running ``geometric-optimize`` and the port number. The worker node must have the QC software installed with the environment variables properly set when ``work_queue_worker`` is executed; one common approach is to write a batch script to execute workers on clusters managed by Slurm or similar job schedulers. If successful, the worker will establish a connection to the master and begin to accept gradient jobs. diff --git a/geometric/neb.py b/geometric/neb.py index 16c2eebf..c84f26af 100644 --- a/geometric/neb.py +++ b/geometric/neb.py @@ -56,7 +56,7 @@ def CoordinateSystem(M, coordtype, chain=False, guessw=0.1): chain : bool True will return a chain object guessw : float - Guessed value for initial Hessians + Guessed weight value for the chain coordinate Returns ------- @@ -965,7 +965,7 @@ def GetCartesianGradient(component, projection): if rmsGrad[n] < self.params.avgg and maxGrad[n] < self.params.maxg: newLocks[n] = True if False not in newLocks: - # HP: In case all of the images are locked before NEB converges, unlock a few. + # HP: In case all the images are locked before NEB converges, unlock a few. logger.info( "All the images got locked, unlocking some images with tighter average gradient value. \n" ) @@ -1635,7 +1635,7 @@ def OptimizeChain(chain, engine, params): # | Adjust Trust Radius / Reject Step |# # =======================================# - chain, trust, trustprint, Y, GW, GP, good = qualitycheck(chain_prev, chain,trust, Quality, ThreLQ, ThreRJ, + chain, trust, trustprint, Y, GW, GP, good = qualitycheck(chain_prev, chain, trust, Quality, ThreLQ, ThreRJ, ThreHQ, Y, GW, GP, Y_prev, GW_prev, GP_prev, params.tmax) if not good: continue diff --git a/geometric/params.py b/geometric/params.py index 11397a62..ff7208cd 100644 --- a/geometric/params.py +++ b/geometric/params.py @@ -477,23 +477,24 @@ def parse_neb_args(*args): grp_univ.add_argument('--nt', type=int, help='Specify number of threads for running in parallel\n(for TeraChem this should be number of GPUs)') grp_nebparam = parser.add_argument_group('nebparam', 'Control the NEB calculation') - grp_nebparam.add_argument('--maxg', type=float, help='Converge when maximum RMS-gradient for any image falls below this threshold (default 0.05 ev/Ang).\n ') - grp_nebparam.add_argument('--avgg', type=float, help='Converge when average RMS-gradient falls below this threshold (default 0.025 ev/Ang).\n ') + grp_nebparam.add_argument('--maxg', type=float, help='Converge when the maximum RMS-gradient of all images falls below this threshold (default 0.05 ev/Ang).\n ') + grp_nebparam.add_argument('--avgg', type=float, help='Converge when the average RMS-gradient of all images falls below this threshold (default 0.025 ev/Ang).\n ') grp_nebparam.add_argument('--guessk', type=float, help='Guess Hessian eigenvalue for displacements (default 0.05).\n ') - grp_nebparam.add_argument('--guessw', type=float, help='Guess weight for chain coordinates (default 0.1).\n ') + #HP 5/10/2024: guessw will be enabled once IC NEB is implemented. + #grp_nebparam.add_argument('--guessw', type=float, help='Guess weight for chain coordinates (default 0.1).\n ') grp_nebparam.add_argument('--nebk', type=float, help='NEB spring constant in units of kcal/mol/Ang^2 (default 1.0).\n ') grp_nebparam.add_argument('--neb_history', type=int, help='Chain history to keep in memory; note chains are very memory intensive, >1 GB each (default 1).\n ') grp_nebparam.add_argument('--neb_maxcyc', type=int, help='Maximum number of chain optimization cycles to perform (default 100).\n ') - grp_nebparam.add_argument('--climb', type=float, help='Activate climbing image for max-energy points when max gradient falls below this threshold (default 0.5).\n ') + grp_nebparam.add_argument('--climb', type=float, help='Activate climbing image for max-energy points when max RMS-gradient falls below this threshold (default 0.5).\n ') grp_nebparam.add_argument('--ncimg', type=int, help='Number of climbing images to expect (default 1).\n ') grp_nebparam.add_argument('--images', type=int, help='Number of NEB images to use (default 11).\n ') grp_nebparam.add_argument('--plain', type=int, help='1: Use plain elastic band for spring force. 2: Use plain elastic band for spring AND potential (default 0).\n ') grp_nebparam.add_argument('--optep', type=str2bool, help='Provide "yes" to optimize two end points of the initial input chain.\n ') grp_nebparam.add_argument('--align', type=str2bool, help='Align images before starting the NEB method (default yes). If "--optep" is set to "yes", the images will be aligned after optimizing the end points.\n ') - grp_nebparam.add_argument('--trust', type=float, help='Starting trust radius, defaults to 0.1 (energy minimization) or 0.01 (TS optimization).\n ') - grp_nebparam.add_argument('--tmax', type=float, help='Maximum trust radius, defaults to 0.3 (energy minimization) or 0.03 (TS optimization).\n ') - grp_nebparam.add_argument('--tmin', type=float, help='Minimum trust radius, do not reject steps trust radius is below this threshold (method-dependent).\n ') - grp_nebparam.add_argument('--skip', action='store_true', help='Skip Hessian updates that would introduce negative eigenvalues.\n ') + grp_nebparam.add_argument('--trust', type=float, help='Starting trust radius (default 0.1). \n ') + grp_nebparam.add_argument('--tmax', type=float, help='Maximum trust radius (default 0.3).\n ') + grp_nebparam.add_argument('--tmin', type=float, help='Minimum trust radius, do not reject steps trust radius is below this threshold.\n ') + grp_nebparam.add_argument('--skip', type=str2bool, help='Skip Hessian updates that would introduce negative eigenvalues.\n ') grp_nebparam.add_argument('--epsilon', type=float, help='Small eigenvalue threshold for resetting Hessian, default 1e-5.\n ') grp_nebparam.add_argument('--port', type=int, help='Work Queue port used to distribute singlepoint calculations. Workers must be started separately.\n')