*********************
Geometry Optimisation
*********************

DL-FIND is a powerful and flexible geometry optimisation library and is the
standard optimiser module in ChemShell. The code is under active development
but the features documented here are those that are available in the version
distributed with ChemShell. 

DL-FIND supports a wide variety of optimisation methods. Standard local
minimisation techniques are supported (steepest descent, conjugate gradient,
Newton-Raphson, BFGS) as well as damped dynamics and limited memory BFGS
(L-BFGS). Transition state methods include nudged elastic band (NEB) for
finding reaction paths and P-RFO and the dimer method for refining transition
states. Conical intersection optimisations are supported using penalty
function, gradient projection and Lagrange-Newton algorithms. Parallel
optimisation methods include a genetic algorithm, stochastic search and
parallel NEB.

DL-FIND accepts input structures in Cartesian coordinates. These may be
converted into a number of internal coordinate systems for the optimisation,
including delocalised internal coordinates (DLCs), total connection, and
hybrid delocalised internal coordinates (HDLCs) for large systems. Constrained
coordinates are supported as well as a general mechanism for imposing
restraints.

A detailed description of DL-FIND's features can be found in:

Johannes Kästner, Joanne M. Carr, Thomas W. Keal, Walter Thiel, Adrian Wander
and Paul Sherwood, J. Phys. Chem. A, 2009, 113 (43), 11856-11865.
DOI: `10.1021/jp9028968 <http://dx.doi.org/10.1021/jp9028968>`_

Please cite this paper if you use DL-FIND in published work.


Running an optimisation
=======================

.. py:class:: Opt

In ChemShell, DL-FIND optimisations can be run using the :py:obj:`Opt()` driver, which
expects a :ref:`Theory <theories>` object, e.g.::

  my_theory = NWChem(...)
  my_job = Opt(theory=my_theory)

The geometry optimisation calculation is run with::

  my_job.run()

Then the results will be stored in :py:obj:`my_job.result` (see: `\<class Result\>`_).

.. _`\<class Result\>`: result.html


Methods
=======

.. py:method:: Opt.run(dryrun=False)

   :param dryrun: If the function is called with ``dryrun=True``, the job will only
                       generate input scripts without executing real calculatons
   :type dryrun: bool, optional

General options
===============

DL-FIND optimisations take the following arguments:

.. py:attribute:: theory
   :noindex:

   Choose the theory level by specifying a :ref:`Theory <theories>` instance

.. py:attribute:: frag2

   Initial second configuration of the system. End point of the initial
   path in NEB calculations, direction of the dimer in dimer calculations.
   Can also be a list of fragments, in which case it defines points along
   the initial NEB path (the desired endpoint should be given last).

.. py:attribute:: maxcycle

   (default: ``100``) Maximum number of optimisation cycles

.. py:attribute:: maxene

   (default: ``10000``) Maximum number of energy evaluations


Coordination system options
===========================

.. py:attribute:: coordinates

   Specifies which coordinate system to use

   Allowed values:

   * ``'cartesian'``: (default) Cartesian coordinates

   * ``'mass'``: Mass-weighted Cartesian coordinates

   * ``'dlc'``: Delocalised internal coordinates (DLC). A redundant set of bonds, angles, torsions, and sometimes inversions is created. A non-redundant combination of them is found by diagonalising the spectroscopic G-matrix. The optimisation is performed in this non-redundant set. `

   * ``'tc'``: Same as ``'dlc'``, but the redundant set consists only of bond - all atoms in the system are connected (total connection scheme)

   * ``'hdlc'``: Hybrid delocalised internal coordinates. The system is partitioned into fragments. Delocalised coordinates (as in DLC) are used within each fragment. The fragments are coupled via Cartesian coordinates. This version is recommended for large systems, as the coordinate transformation scales linearly with the number of atoms. 

   * ``'hdlc-tc'``: Hybrid delocalised internal coordinates based on the total connection scheme. 

:NB: that the ``'dlc'`` option is not equivalent to specifying hdlc with a 
   single residue. The ``'dlc'`` case gives pure internal delocalised coordinates 
   (6 external degrees of freedom are dropped), whereas with ``'hdlc'`` Cartesian 
   coordinates are added to the primitive list even if there is only one 
   residue, and the external degrees of freedom are retained.

.. py:attribute:: residues

   Defines residues (atom groups) for use with HDLC coordinates.

   Accepts a dictionary of the form::

     { 'res1' : [atom1, atom2, atom3, ...],
       'res2' : [atom4, atom5, atom6, ...] }

   NB: For biomolecular Fragments this dictionary can be generated from PDB data using
   the getHDLCResidues() method on the fragment.

.. py:attribute:: constraints

   Define constraints for use with either HDLC or Cartesian coordinates.

   For HDLC coordinates, internal constraints can be defined using a list of the form::

     [ ('bond', atom1, atom2),
       ('angle', atom1, atom2, atom3),
       ('torsion', atom1, atom2, atom3, atom4),
      ... ]

   For Cartesian coordinates, constraints can be defined using a list of the form::

     [ ('x', atom1), ('yz', atom2), ... ]

   where the first constraint freezes the X coordinate of atom1, the second freezes 
   the Y and Z coordinates of atom2, etc.
 
  
Optimisation algorithm options
==============================

.. py:attribute:: algorithm

   Specifies the optimisation algorithm

   Allowed values:
    * ``'lbfgs'``: (default) Limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) optimisation.
      Recommended for minimum search, NEB, and the dimer method. The time and
      memory requirements spent for determining the search direction and step 
      length scale linearly with the system size. The number of steps kept in 
      memory is specified by lbfgs_mem.
    * ``'prfo'``: The partitioned rational function optimisation method. 
      Recommended for transition state search (unless done by the dimer or NEB 
      methods). It requires the calculation of the Hessian. The keyword 
      update_method can be used to specify how the Hessian should be updated.
    * ``'cg'``: Conjugate gradient method following Polak–Ribiere.
    * ``'sd'``: Steepest descent.
    * ``'dyn'``: Damped molecular dynamics. Four additional parameters can be 
      specified. timestep is the time step in atomic units. :py:attr:`fric0` is the 
      start friction, :py:attr:`fricfac` is the factor to reduce the friction each 
      time the energy decreases, and :py:attr:`fricp` is the friction to apply if the 
      energy increases. The frictions are defined so that ``0`` corresponds to 
      free (undamped) dynamics, and ``1`` corresponds to steepest descent. 
      This variable friction facilitates convergence to an energy minimum.
    * ``'nr'``: Basic Newton-Raphson/quasi-Newton optimiser.
    * ``'ln'``: Lagrange-Newton constrained optimiser (for conical 
      intersection searches only). 

.. py:attribute:: trust_radius

   Specifies the trust radius approach

   Allowed values:

   * ``'constant'``: (default) Constant trust radius (equal to maxstep)

   * ``'energy'``: Trust radius based on the energy change

   * ``'gradient'``: Trust radius based on the projection of the gradient on the current step. Can be used when no energy expression is available (e.g. NEB).

.. py:attribute:: maxstep

   (Default: ``0.5``) Maximum change of one coordinate component in one 
   optimisation step.

.. py:attribute:: tolerance

   (Default: ``0.00045``) Convergence criterion for the maximum gradient 
   component. The criterion for the maximum step component is :py:attr:`4*tolerance`,
   for the RMS step it is :py:attr:`8/3*tolerance`, for the RMS gradient it is
   :py:attr:`2/3*tolerance`. All criteria are applied in the coordinate system 
   specified by ``coordinates``.

.. py:attribute:: tolerance_e

   (Default: ``tolerance/450``) Convergence criterion for the energy 
   change (in Hartree).


Nudged elastic band options
===========================

.. py:attribute:: neb

   Specifies if and how the nudged elastic band (NEB) method is to be used

   Allowed values:

   * ``'no'``: (default) No NEB calculation will be performed. :py:attr:`neb` is automatically set to ``'no'`` if no :py:attr:`frag2` is specifed.

   * ``'free'``: The endpoints of the NEB path are minimised

   * ``'frozen'``: The endpoints of the NEB path are frozen

   * ``'perpendicular'``: The endpoints of the NEB path minimised perpendicular to the NEB path. Their distance to the nearest other point is kept fixed. 

.. py:attribute:: nimages

   (Default: ``1``) Specifies the number of images in a NEB calculation.

.. py:attribute:: nebk

   (Default: ``0.01``) Force constant for NEB calculations.


Dimer method options
====================

.. py:attribute:: dimer

   Specifies if the dimer method for transition state (TS) search is 
   to be used.

   Allowed values:

   * ``False``: (default) Do not use the dimer method
   * ``True``: Use the dimer method for TS search

.. py:attribute:: delta

   (Default: ``0.01``) Dimer separation (in internal coordinates)

.. py:attribute:: tsrelative

   (Default: ``False``) If True, the coordinates given in frag2 are given as displacements relative to the coordinates of frag (i.e. the frag2 coords will
   equal frag.coords + frag2.coords).


.. TODO lots more options etc. in above tables and others
.. TODO so far I have only included options used in the runs in tests/