
.. include:: autodoc_abbr_options_c.rst

.. index::
   single: geometry optimization
   see: optimization; geometry optimization

.. _`sec:optking`:

Geometry Optimization
=====================

.. codeauthor:: Rollin A. King
.. sectionauthor:: Rollin A. King and Lori A. Burns

*Module:* :ref:`Keywords <apdx:optking>`, :ref:`PSI Variables <apdx:optking_psivar>`, :source:`OPTKING <src/bin/optking>`

|PSIfour| carries out molecular optimizations using a module called
optking.  The optking program takes as input nuclear gradients and,
optionally, nuclear second derivatives |w---w| both in Cartesian coordinates.
The default minimization algorithm employs an empirical model Hessian,
redundant internal coordinates, a RFO step, and the BFGS Hessian update.

The principal literature references include the introduction of redundant
internal coordinates by Peng et al. [Peng:1996:49]_.
The general approach employed in this code
is similar to the "model Hessian plus RF method" described and tested by Bakken and
Helgaker [Bakken:2002:9160]_. (However, for separated
fragments, we have chosen not to employ by default their "extra-redundant"
coordinates defined by their "auxiliary interfragment" bonds.  These can be
included via the option |optking__add_auxiliary_bonds|).

The internal coordinates are generated automatically based on an assumed bond
connectivity.  The connectivity is determined by testing if the interatomic
distance is less than the sum of atomic radii times the value of
|optking__covalent_connect|. If the user finds that some
connectivity is lacking by default, then this value may be increased.
Otherwise, the internal coordinate definitions may be modified.  If one
desires to see or modify the internal coordinates being used, then one can set
|optking__intcos_generate_exit| to true.  The internal coordinate
definitions are provided in the file with extension ".intco".  See the :ref:`sec:optkingExamples`
section for more detail.

.. warning:: Optimizations where the molecule is specified in Z-matrix format 
   with dummy atoms will result in the molecule being converted to a Cartesian representation.

The ongoing development of optking is providing for unique treatment of
coordinates which connect distinct molecular fragments.  Thus, several keywords
relate to "interfragment modes", though many of these capabilities are
still under development.  Presently by default, separate fragments are bonded by
nearest atoms, and the whole system is treated as if it were part of one
molecule.  However, with the option |optking__frag_mode|, fragments
may instead be related by a unique set of interfragment coordinates defined by
reference points within each fragment.  The reference points can be atomic
positions (current default), linear combinations of
atomic positions, or located on the principal axes (not yet working).

Basic Keywords
^^^^^^^^^^^^^^

.. include:: autodir_options_c/optking__opt_type.rst
.. include:: autodir_options_c/optking__step_type.rst
.. include:: autodir_options_c/optking__geom_maxiter.rst
.. include:: autodir_options_c/optking__g_convergence.rst
.. include:: autodir_options_c/optking__full_hess_every.rst
.. include:: autodir_options_c/optking__intcos_generate_exit.rst

.. index:: geometry optimization; minima
.. _`sec:optkingExamples`:

Optimizing Minima
^^^^^^^^^^^^^^^^^

First, define the molecule and basis in the input. ::

   molecule h2o {
     O
     H 1 1.0
     H 1 1.0 2 105.0
   }
    
   set basis dz

Then the following are examples of various types of calculations that can be completed.

* Optimize a geometry using default methods (RFO step)::

   optimize('scf')

* Optimize using Newton-Raphson steps instead of RFO steps::

   set step_type nr
   optimize('scf')

* Optimize using energy points instead of gradients::

   optimize('scf', dertype='energy')

* Optimize while limiting the initial step size to 0.1 au::

   set intrafrag_step_limit 0.1
   optimize('scf')

* Optimize while always limiting the step size to 0.1 au::

   set {
     intrafrag_step_limit     0.1
     intrafrag_step_limit_min 0.1
     intrafrag_step_limit_max 0.1
   }

   optimize('scf')

* Optimize while calculating the Hessian at every step::

   set full_hess_every 1
   optimize('scf')

Hessian
^^^^^^^

If Cartesian second derivatives are available, optking can read them
and transform them into internal coordinates to make an initial Hessian in
internal coordinates.  Otherwise, several empirical Hessians are available,
including those of Schlegel [Schlegel:1984:333]_ and Fischer and Almlof
[Fischer:1992:9770]_.
Either of these or a simple diagonal Hessian may be selected using the 
|optking__intrafrag_hess| keyword.

All the common Hessian update schemes are available.  For formulas, see
Schlegel [Schlegel:1987:AIMQC]_ and Bofill [Bofill:1994:1]_.

The Hessian may be computed during an optimization using the 
|optking__full_hess_every| keyword.

.. index:: 
   pair: geometry optimization; transition state
   pair: geometry optimization; IRC
   single: geometry optimization; constrained

Transition States, Reaction Paths, and Constrained Optimizations
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

* Calculate a starting Hessian and optimize the "transition state" of
  linear water (note that without a reasonable starting geometry and
  Hessian, such a straightforward search often fails)::

   molecule h2o {
      O
      H 1 1.0
      H 1 1.0 2 160.0
   }

   set {
      basis dz
      full_hess_every 0
      opt_type ts
   }

   optimize('scf')

* At a transition state (planar HOOH), compute the second derivative, and
  then follow the intrinsic reaction path to the minimum::

   molecule hooh {
      symmetry c1
      H
      O 1 0.946347
      O 2 1.397780 1  107.243777
      H 3 0.946347 2  107.243777   1 0.0
   }

   set {
      basis dzp
      opt_type irc
      geom_maxiter 50
   }

   frequencies('scf')
   optimize('scf')

* Generate the internal coordinates and then stop::

   set intcos_generate_exit true
   optimize('scf')

  The coordinates may then be found in the "intco" file.  In this case, the file contains::
  
     F 1 3
     R      1     2
     R      1     3
     B      2     1     3
  
  The first line indicates a fragment containing atoms 1-3.  The following lines define
  two distance coordinates (bonds) and one bend coordinate.  This file can be modified, and if present,
  is used in subsequent optimizations.  Since the multiple-fragment coordinates are still under
  development, they are not documented here.  However, if desired, one can change the value
  of |optking__frag_mode|, generate the internal coordinates, and see how multiple
  fragment systems are defined.
  
  Coordinates may be frozen or fixed by adding an asterisk after the letter of the coordinate.
  To optimize with the bond lengths fixed at their initial values, it is currently necessary to
  generate and then modify the internal coordinate definitions as follows::
  
     F 1 3
     R*     1     2
     R*     1     3
     B      2     1     3

.. index:: geometry optimization; convergence criteria

Convergence Criteria
^^^^^^^^^^^^^^^^^^^^

Optking monitors five quantities to evaluate the progress of a geometry 
optimization. These are (with their keywords) the change in energy 
(|optking__max_energy_g_convergence|), the maximum element of 
the gradient (|optking__max_force_g_convergence|), the root-mean-square 
of the gradient (|optking__rms_force_g_convergence|), the maximum element
of displacement (|optking__max_disp_g_convergence|), and the 
root-mean-square of displacement (|optking__rms_disp_g_convergence|), 
all in internal coordinates and atomic units. Usually, these options will not 
be set directly. Primary control for geometry convergence lies with the keyword 
|optking__g_convergence| which sets the aforementioned in accordance 
with Table :ref:`Geometry Convergence <table:optkingconv>`.

|
|

.. _`table:optkingconv`:

.. table:: Summary of sets of geometry optimization criteria available in |PSIfour|

    +-----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+
    | |optking__g_convergence|    | Max Energy                 | Max Force                  | RMS Force                  | Max Disp                   | RMS Disp                   |
    +=============================+============================+============================+============================+============================+============================+
    | NWCHEM_LOOSE [#fd]_         |                            | :math:`4.5 \times 10^{-3}` | :math:`3.0 \times 10^{-3}` | :math:`5.4 \times 10^{-3}` | :math:`3.6 \times 10^{-3}` |
    +-----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+
    | GAU_LOOSE [#ff]_            |                            | :math:`2.5 \times 10^{-3}` | :math:`1.7 \times 10^{-3}` | :math:`1.0 \times 10^{-2}` | :math:`6.7 \times 10^{-3}` |
    +-----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+
    | TURBOMOLE [#fd]_            | :math:`1.0 \times 10^{-6}` | :math:`1.0 \times 10^{-3}` | :math:`5.0 \times 10^{-4}` | :math:`1.0 \times 10^{-3}` | :math:`5.0 \times 10^{-4}` |
    +-----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+
    | GAU [#fc]_ [#ff]_           |                            | :math:`4.5 \times 10^{-4}` | :math:`3.0 \times 10^{-4}` | :math:`1.8 \times 10^{-3}` | :math:`1.2 \times 10^{-3}` |
    +-----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+
    | CFOUR [#fd]_                |                            |                            | :math:`1.0 \times 10^{-4}` |                            |                            |
    +-----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+
    | QCHEM [#fa]_ [#fe]_         | :math:`1.0 \times 10^{-6}` | :math:`3.0 \times 10^{-4}` |                            | :math:`1.2 \times 10^{-3}` |                            |
    +-----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+
    | MOLPRO [#fb]_ [#fe]_        | :math:`1.0 \times 10^{-6}` | :math:`3.0 \times 10^{-4}` |                            | :math:`3.0 \times 10^{-4}` |                            |
    +-----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+
    | GAU_TIGHT [#fc]_ [#ff]_     |                            | :math:`1.5 \times 10^{-5}` | :math:`1.0 \times 10^{-5}` | :math:`6.0 \times 10^{-5}` | :math:`4.0 \times 10^{-5}` |
    +-----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+
    | GAU_VERYTIGHT [#ff]_        |                            | :math:`2.0 \times 10^{-6}` | :math:`1.0 \times 10^{-6}` | :math:`6.0 \times 10^{-6}` | :math:`4.0 \times 10^{-6}` | 
    +-----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+----------------------------+

.. rubric:: Footnotes

.. [#fa] Default
.. [#fb] Baker convergence criteria are the same.
.. [#fc] Counterpart NWCHEM convergence criteria are the same.
.. [#fd] Convergence achieved when all active criteria are fulfilled.
.. [#fe] Convergence achieved when **Max Force** and one of **Max Energy** or **Max Disp** are fulfilled.
.. [#ff] Normal convergence achieved when all four criteria (**Max Force**, **RMS Force**,
         **Max Disp**, and **RMS Disp**) are fulfilled. To help with flat 
         potential surfaces, alternate convergence achieved when 100\ :math:`\times`\ *rms force* is less 
         than **RMS Force** criterion.

For ultimate control, specifying a value for any of the five monitored options activates that
criterium and overwrites/appends it to the criteria set by |optking__g_convergence|.
Note that this revokes the special convergence arrangements detailed in notes [#fe]_ and [#ff]_ 
and instead requires all active criteria to be fulfilled to 
achieve convergence. To avoid this revokation, turn on keyword |optking__flexible_g_convergence|.

.. index::
   pair: geometry optimization; output

Output
^^^^^^

The progress of a geometry optimization can be monitored by grepping the output file for the
tilde character (``~``). This produces a table like the one below that shows
for each iteration the value for each of the five quantities and whether the criterion
is active and fulfilled (``*``), active and unfulfilled ( ),  or inactive (``o``). ::

   --------------------------------------------------------------------------------------------- ~
    Step     Total Energy     Delta E     MAX Force     RMS Force      MAX Disp      RMS Disp    ~
   --------------------------------------------------------------------------------------------- ~
     Convergence Criteria    1.00e-06 *    3.00e-04 *             o    1.20e-03 *             o  ~
   --------------------------------------------------------------------------------------------- ~
       1     -38.91591820   -3.89e+01      6.91e-02      5.72e-02 o    1.42e-01      1.19e-01 o  ~
       2     -38.92529543   -9.38e-03      6.21e-03      3.91e-03 o    2.00e-02      1.18e-02 o  ~
       3     -38.92540669   -1.11e-04      4.04e-03      2.46e-03 o    3.63e-02      2.12e-02 o  ~
       4     -38.92548668   -8.00e-05      2.30e-04 *    1.92e-04 o    1.99e-03      1.17e-03 o  ~
       5     -38.92548698   -2.98e-07 *    3.95e-05 *    3.35e-05 o    1.37e-04 *    1.05e-04 o  ~

The full list of keywords for optking is provided in Appendix :ref:`apdx:optking`.

Information on the Psithon function that drives geometry optimizations is provided
at :py:func:`~driver.optimize`.

