<set title="DL-FIND">
<set tag="opt">
<set tag1="dlfind">

<file include="head.htp">

<h2><a name="intro" />1. Introduction</h2>

<p>DL-FIND is a powerful and flexible geometry optimisation library and 
is the recommended 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 can also be downloaded as a standalone library from the <a
href="http://ccpforge.cse.rl.ac.uk/projects/dl-find/">CCPForge repository</a>.</p>

<p>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.</p>

<p>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 <a href="hdlc.html">hybrid delocalised 
internal coordinates</a> (HDLCs) for large systems. Constrained coordinates are supported 
as well as a general mechanism for imposing restraints.</p>

<p>A detailed description of DL-FIND's features can be found in:</p>

<p>Johannes K&auml;stner, Joanne M. Carr, Thomas W. Keal, Walter Thiel, Adrian Wander and Paul Sherwood, <i>J. Phys. Chem. A</i>, 2009, 113 (43), 11856-11865. DOI: <a href="http://dx.doi.org/10.1021/jp9028968">10.1021/jp9028968</a></p>

<p>Please cite this paper if you use DL-FIND in published work.</p>

<h2><a name="cmd" />2. Command line arguments</h2>

<h3><a name="general" />2.1 General options</h3>

<TABLE CELLPADDING=3 BORDER="1">
<TR><TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw1><I>Argument</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw2><I>Argument type</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw3><I>Mandatory</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw4><I>Default</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw5><I>To specify</I></TD>
</TR>

<TR>
<TD ALIGN="LEFT" VALIGN="TOP">theory</TD>
<TD ALIGN="LEFT" VALIGN="TOP">string</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">mndo</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Module used to compute energy and forces</TD>
</TR>

<TR>
<TD ALIGN="LEFT" VALIGN="TOP">coords</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use fragment_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">yes</TD>
<TD ALIGN="LEFT" VALIGN="TOP">none </TD>
<TD ALIGN="LEFT" VALIGN="TOP">Initial configuration of the system</TD>
</TR>

<TR>
<TD ALIGN="LEFT" VALIGN="TOP">result</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use fragment_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">dl-find.result</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Optimised coordinates</TD>
</TR>

<TR>
<TD ALIGN="LEFT" VALIGN="TOP">coords2</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use fragment_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">yes (for NEB)</TD>
<TD ALIGN="LEFT" VALIGN="TOP">none </TD>
<TD ALIGN="LEFT" VALIGN="TOP">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. Then it defines points along the
initial NEB path (the desired endpoint should be given last).</TD>
</TR>

<TR>
<TD ALIGN="LEFT" VALIGN="TOP">result2</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use fragment_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">none</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Optimised second coordinate set. Its contents
depend on the task. Usually a transition mode.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">tsrelative</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">false</TD>
<TD ALIGN="LEFT" VALIGN="TOP">coords2 and result2 are specified relative to
coords and result, respectively</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">maxcycle</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">100</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Maximum number of optimisation cycles</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">maxene</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">10000</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Maximum number of energy and gradient evaluations</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">dump</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Write restart information every <tt>dump</tt>
steps (0: never).</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">restart</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">false</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Restart from a dump file.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">list_option</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use listopt_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">medium</TD>
<TD ALIGN="LEFT" VALIGN="TOP">How much output to generate (including generation of .xyz files)</TD>
</TR>
</TABLE>

<p>Setting <tt>list_option=full</tt> for a standard minimisation will output a 
<tt>path.xyz</tt> file containing the optimisation trajectory and 
<tt>path_active.xyz</tt> containing only the active atoms. Other methods 
will generate other .xyz files as appropriate (see below).</p>

<h3><a name="coord" />2.2 Coordinate system options</h3>

<TABLE CELLPADDING=3 BORDER="1">
<TR><TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw1><I>Argument</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw2><I>Argument type</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw3><I>Mandatory</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw4><I>Default</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw5><I>To specify</I></TD>
</TR>

<TR>
<TD ALIGN="LEFT" VALIGN="TOP">coordinates</TD>
<TD ALIGN="LEFT" VALIGN="TOP">keyword</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">cartesian</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Which coordinate system to use. Can be
<tt>cartesian</tt>,  <tt>mass</tt>,  <tt>dlc</tt>,  <tt>tc</tt>,
<tt>hdlc</tt>, or  <tt>hdlc-tc</tt>.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">residues</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use list_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">none</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Definition of residues (fragments) for the option
<tt>coordinates=hdlc</tt>. <A HREF="hdlc.html#residues">See HDLCs</a>.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">frozen</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use list_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">none</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Specification of the atoms to be
frozen. Complementary to <tt>active_atoms</tt> </TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">active_atoms</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use list_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">all atoms</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Specification of the atoms to be optimised.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">constraints</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use list_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">none</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><A HREF="hdlc.html#constraints">Constraints</A> specification for use with internal coordinates.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">restraints</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use list_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">none</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><a href="tools.html#restraints">Restraints</a> to apply</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">connect</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use list_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">none</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Add user-defined connections, format: { atom1
atom2 }</TD>
</TR>
</TABLE>

<p>The keyword <tt>coordinates</tt> specifies the coordinate system in which the
optimisation should be performed. Possible choices are:</p>
<ul>
<li> <tt>cartesian</tt>: Cartesian coordinates 
<li> <tt>mass</tt>: Mass-weighted Cartesian coordinates
<li> <tt>dlc</tt>: Delocalised internal coordinates. 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.
<li> <tt>tc</tt>: Same as <tt>dlc</tt>, but the redundant set consists only of bonds&mdash;all
  atoms in the system are connected (total connection scheme).
<li> <tt>hdlc</tt>: <a href="hdlc.html">Hybrid delocalised internal coordinates</a>. 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.
<li> <tt>hdlc-tc</tt>:  Hybrid delocalised internal coordinates based on the
  total connection scheme.
</ul>

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

<h3><a name="opt" />2.3 Optimisation algorithm options</h3>

<TABLE CELLPADDING=3 BORDER="1">
<TR><TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw1><I>Argument</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw2><I>Argument type</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw3><I>Mandatory</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw4><I>Default</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw5><I>To specify</I></TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">optimiser</TD>
<TD ALIGN="LEFT" VALIGN="TOP">keyword</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">lbfgs</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Specifies the optimisation algorithm. Can be
<tt>lbfgs</tt>, <tt>prfo</tt>, <tt>cg</tt>, <tt>sd</tt>,
<tt>dyn</tt>, <tt>nr</tt> or <tt>ln</tt>.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">trustradius</TD>
<TD ALIGN="LEFT" VALIGN="TOP">keyword</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">constant</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Specifies the trust radius approach. Can be
<tt>constant</tt>, <tt>energy</tt>, or <tt>gradient</tt>.</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">maxstep</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.5</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Maximum change of one coordinate component in one
optimisation step.</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">tolerance</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.00045</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Convergence criterion for the maximum gradient
component. The criterion for the maximum step component is 4*<tt>tolerance</tt>,
for the RMS step it is 8/3*<tt>tolerance</tt>, for the RMS gradient it is
2/3*<tt>tolerance</tt>. All criteria are applied in the coordinate system
specified by <tt>coordinates</tt>.</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">tolerance_e</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">tolerance/450</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Convergence criterion for the energy change (in Hartree).</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">lbfgs_mem</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">DOF</TD>
<TD ALIGN="LEFT" VALIGN="TOP">The number of steps kept in the L-BFGS
memory. Default: number of degrees of freedom.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">initial_hessian</TD>
<TD ALIGN="LEFT" VALIGN="TOP">keyword</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">external</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Type of Hessian to calculate at the start or when resetting. Can be <tt>external</tt>, <tt>onepoint</tt>, <tt>twopoint</tt>, <tt>diagonal</tt>, or <tt>identity</tt>.</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">update_method</TD>
<TD ALIGN="LEFT" VALIGN="TOP">keyword</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">bofill</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Update algorithm for the Hessian 
(e.g. in P-RFO TS search or quasi-Newton optimisation). Can be
<tt>bofill</tt>, <tt>powell</tt>, <tt>bfgs</tt>, or <tt>none</tt>.</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">delta</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.01</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Step size in the finite-difference Hessian calculation.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">maxupdate</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">50</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Maximum number of Hessian updates before the
Hessian is recalculated.</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">soft</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.005</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Eigenmodes with an absolute Eigenvalue below
<tt>soft</tt> are ignored by the P-RFO optimiser. This avoids movements in the
translation and rotation direction. A maximum of 6 modes are ignored. No modes
are considered soft when internal coordinates are used.</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">minstep</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.00001</TD>
<TD ALIGN="LEFT" VALIGN="TOP">If the step size is below <tt>minstep</tt>, 
the Hessian is not updated (to avoid problems with noisy gradients).</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">scalestep</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">1.0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Scale the step obtained from SD or CG algorithms.</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">distort</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Distort the starting geometry along coords2 (if
positive) or against it (if negative). This may be useful when starting from a
TS and searching for minima.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">timestep</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">1.0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Time step (in a.u.) for damped molecular dynamics
(<a href="#dyn"><tt>optimiser=dyn</tt></a>).</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">fric0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.3</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Start friction for damped molecular dynamics
(<a href="#dyn"><tt>optimiser=dyn</tt></a>).</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">fricfac</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.95</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Factor to reduce the friction if the energy
decreases in damped molecular dynamics
(<a href="#dyn"><tt>optimiser=dyn</tt></a>).</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">fricp </TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.3</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Friction to apply if the energy increases in
damped molecular dynamics 
(<a href="#dyn"><tt>optimiser=dyn</tt></a>).</TD>  
</TR>
</TABLE>

<p>The keyword <tt>optimiser</tt> specifies the optimisation algorithm. Possible
choices are:</p>

<ul>
<li> <tt>lbfgs</tt>: Limited-memory
  Broyden&ndash;Fletcher&ndash;Goldfarb&ndash;Shanno 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 <tt>lbfgs_mem</tt>.
<li> <tt>prfo</tt>: 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 <tt>update_method</tt> can be used to specify how the
  Hessian should be updated.
<li> <tt>cg</tt>: Conjugate gradient method following Polak&ndash;Ribiere.
<li> <tt>sd</tt>: Steepest descent.
<li> <tt>dyn</tt> <a name="dyn"></a>: Damped molecular dynamics. Four
  additional parameters can be 
  specified. <tt>timestep</tt> is the time step in atomic units. <tt>fric0</tt>
  is the start friction, <tt>fricfac</tt> is the factor to reduce the friction
  each time the energy decreases, and <tt>fricp</tt> 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.
<li> <tt>nr</tt>: Basic Newton-Raphson/quasi-Newton optimiser.
<li> <tt>ln</tt>: Lagrange-Newton constrained optimiser (for conical intersection searches only).
</ul>

<p>The trust radius approach (or line search) is controlled via the keyword
<tt>trustradius</tt>. Possible choices are:</p>

<ul>
<li> <tt>constant</tt>: Constant trust radius (equal to <tt>maxstep</tt>)
<li> <tt>energy</tt>: Trust radius based on the energy change
<li> <tt>gradient</tt>: 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).
</ul>

<p>For optimisers requiring a Hessian, such as the Newton-Raphson or P-RFO algorithms, 
the type of Hessian to be calculated is specified by the keyword
<tt>initial_hessian</tt>. Possible choices are:</p>

<ul>
<li> <tt>external</tt>: Calculate with the ChemShell routine <tt>$theory.hessian</tt> (default)
<li> <tt>onepoint</tt>: Build by one point finite difference of the gradient
<li> <tt>twopoint</tt>: Build by two point finite difference of the gradient 
<li> <tt>diagonal</tt>: Build a diagonal Hessian by a single one point finite difference
<li> <tt>identity</tt>: Set the Hessian to be an identity matrix
</ul>

<p>In the case of <tt>initial_hessian=external</tt>, two point finite difference will be 
used if an external Hessian is not available.</p>

<p>The Hessian will be updated depending on the
keyword <tt>update_method</tt>. Possible choices for <tt>update_method</tt>
are:</p>

<ul>
<li> <tt>bofill</tt>: Bofill update of the Hessian (default)
<li> <tt>powell</tt>: Powell update of the Hessian 
<li> <tt>bfgs</tt>: BFGS update of the Hessian. (Note: this is a full BFGS Hessian update, not to be confused with the limited-memory version).
<li> <tt>none</tt>: No update, recalculate the Hessian in each step
</ul>

<h3><a name="neb" />2.4 Nudged elastic band options</h3>

<TABLE CELLPADDING=3 BORDER="1">
<TR><TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw1><I>Argument</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw2><I>Argument type</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw3><I>Mandatory</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw4><I>Default</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw5><I>To specify</I></TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">neb</TD>
<TD ALIGN="LEFT" VALIGN="TOP">keyword</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Specifies if and how the nudged elastic band method is to
be used. Can be <tt>no</tt>, <tt>free</tt>, <tt>frozen</tt>, or
<tt>perpendicular</tt>.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">nimage</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">1</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Specifies the number of images in a NEB
calculation.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">neb_cart</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">false</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Find the NEB starting path in the coordinate
system specified by coordinates, but switch to Cartesian coordinates for the
optimisation.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">nebk</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.01</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Force constant for NEB calculations.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">neb_climb_test</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">3.0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Threshold for spawning the climbing image (as a 
multiple of the gradient tolerance for convergence).</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">neb_freeze_test</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">1.0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Threshold for freezing images during the optimisation (as a 
multiple of the gradient tolerance for convergence).</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">weights</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use list_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">all 1</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Weights specified as a Tcl list of real
numbers. One number for each atom. 
The higher the weight for an atom, the more it contributes to the NEB path. 
Atoms with weight=0 are excluded from the
transition mode.</TD>
</TR>
</TABLE>

<p>The keyword <tt>neb</tt> specifies that a nudged elastic band calculation
should be performed. The improved-tangent NEB algorithm with a climbing image
is implemented. Three versions are available:</p>

<ul> 
<li> <tt>free</tt>: The endpoints of the NEB path are minimised
<li> <tt>frozen</tt>: The endpoints of the NEB path are frozen
<li> <tt>perpendicular</tt>: The endpoints of the NEB path minimised 
perpendicular to the NEB path. Their distance to the nearest other point is
kept fix.
<li> <tt>neb= no</tt> no NEB is used.
</ul>

<p>The NEB path will be initiated as a linear transit between the geometries
specified by the keywords <tt>coords</tt> and <tt>coords2</tt>. The linear
transit will be calculated in the coordinate system specified by
<tt>coordinates</tt>. Sometimes it is preferable to initiate the path in
internal coordinates (in particular <tt>coordinates=tc</tt>) but to optimise
the path in Cartesian coordinates (more stable). In this case, specify
<tt>neb_cart=true</tt>.</p>

<p>The number of images is specified by <tt>nimage</tt>. The calculation will
start out without climbing image, thus with one image less than
<tt>nimage</tt>. Once the gradient perpendicular to the path is small enough (by default 3
times the general tolerance for gradient components), the climbing image will
be placed at the interpolated maximum of the path. From this point on, it will
be optimised with the rest of the path. Once the perpendicular force of an
image is smaller than the general tolerance for gradient components,
the respective image is frozen. This saves a significant number of energy
evaluations. The thresholds for 
spawning the climbing image and freezing images can be set by <tt>neb_climb_test</tt>
and <tt>neb_freeze_test</tt> respectively. Setting them to 0.0 will switch off these 
features, which may help if you experience convergence problems during the later 
stage of the optimisation.</p>

<p>The use of the L-BFGS optimiser is recommended with NEB.</p>

<p>Additional output files may be written when using NEB. With
<tt>list_option=medium</tt> a file called nebinfo will be written after each
completed NEB step. It contains the path length to each image, the energy of
each image relative to the first image, and the work (gradient projected on the
path vector). The latter may be more relevant than the energy for huge
systems. Additionally a file nebpath.xyz is written. It contains the
present geometry of each image. ChemShell fragment files neb_*.c with the same
information are also written. With <tt>list_option=full</tt> additional files
neb_*.xyz are written containing the optimisation trajectory for each
image. These can be viewed e.g. with "vmd -m neb_*.xyz".</p>

<p>The NEB method may be run in parallel 
using the <a href="parallel.html#taskfarm">task-farming framework</a>.</p>

<h3><a name="dimer" />2.5 Dimer method options</h3>

<TABLE CELLPADDING=3 BORDER="1">
<TR><TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw1><I>Argument</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw2><I>Argument type</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw3><I>Mandatory</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw4><I>Default</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw5><I>To specify</I></TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">dimer</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">false</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Specifies if the dimer method for TS search is to
be used.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">dimer_interpolate</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">true</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Specifies if the gradient in the dimer rotation
should be interpolated.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">maxrot</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">10</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Maximum number of dimer rotation steps before a
translation step.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">tolrot</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">5.0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Tolerance for dimer rotation (in degrees).</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">delta</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.01</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Dimer separation (in internal coordinates).</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">weights</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use list_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">all 1</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Weights specified as a Tcl list of real
numbers. One number for each atom. The higher the weight for an atom, 
the more it contributes to the
transition mode. Atoms with weight=0 are excluded from the
transition mode.</TD>
</TR>
</TABLE>

<p>The keyword <tt>dimer</tt> is used to start a transition state search using the dimer
method. Two images of the system (their distance is
specified by <tt>delta</tt>) are calculated. They are optimised along the
force in all directions perpendicular to the dimer axis, and against the force
along the dimer axis. Thus the system converges to a first order saddle point.</p>

<p>The translation of the dimer is optimised with the algorithm specified by
<tt>optimiser</tt> (L-BFGS is recommended). The rotation is optimised by a
conjugate gradient method. 
After a rotation step, the gradient can be interpolated (default) or
recalculated (<tt>dimer_interpolate=false</tt>).</p>

<h3><a name="conint" />2.6 Conical intersection options</h3>

<TABLE CELLPADDING=3 BORDER="1">
<TR><TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw1><I>Argument</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw2><I>Argument type</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw3><I>Mandatory</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw4><I>Default</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw5><I>To specify</I></TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">conint</TD>
<TD ALIGN="LEFT" VALIGN="TOP">keyword</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Specifies if and how multi-state calculations are carried out.
Can be <tt>no</tt>, <tt>pf</tt>, <tt>gp</tt>, or
<tt>ln</tt>.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">state_i</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">1</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Choice of lower state for conical intersection optimisation.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">state_j</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">2</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Choice of upper state for conical intersection optimisation.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">coupled_states</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">yes</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Specify whether an interstate coupling gradient exists (if false, the coupling is set to zero). Only applicable to the <tt>gp</tt> and <tt>ln</tt> methods.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">pf_c1</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">5.0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">First weighting option for the penalty function term in the penalty 
function conical intersection optimisation algorithm (in kcal/mol).</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">pf_c2</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">5.0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Second weighting option for the penalty function term in the penalty 
function conical intersection optimisation algorithm (in kcal/mol).</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">gp_c3</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">1.0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Overall scaling of the gradient in the gradient projection
conical intersection optimisation algorithm.</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">gp_c4</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.9</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Relative weighting of the two terms of the gradient in the 
gradient projection conical intersection optimisation algorithm.</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">ln_t1</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.0001</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Threshold for the Lagrange-Newton conical intersection 
optimisation algorithm (in kcal/mol), below which orthogonalisation is switched on.</TD>  
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">ln_t2</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">1.0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Threshold for the Lagrange-Newton conical intersection 
optimisation algorithm (in kcal/mol), above which orthogonalisation is switched off.</TD>  
</TR>
</table>

<p>The keyword <tt>conint</tt> specifies that a conical intersection
optimisation should be carried out. At each cycle energies, gradients, and (if
applicable) interstate coupling gradients are calculated for two electronic
states. This information is used to search for a local minimum energy crossing
point between the states according to the chosen algorithm. The algorithms
available are:</p> 

<ul> 
<li> <tt>pf</tt>: Penalty function method
<li> <tt>gp</tt>: Gradient projection method
<li> <tt>ln</tt>: Lagrange-Newton method
<li> <tt>no</tt>: no multi-state calculation (default).
</ul>

<p>The states to be optimised are chosen using the options <tt>state_i</tt> and
<tt>state_j</tt>. Specific algorithm options can be set using <tt>pf_c1</tt>
and <tt>pf_c2</tt> for the penalty function method, <tt>gp_c3</tt> and
<tt>gp_c4</tt> for the gradient projection method, and <tt>ln_t1</tt> and 
<tt>ln_t2</tt> for the Lagrange-Newton method. The default values are those
recommended in the <a href="http://dx.doi.org/10.1007/s00214-007-0331-5">MNDO
implementation of these algorithms</a>.</p> 

<p>By default the gradient projection and Lagrange-Newton methods require the calling program to provide the interstate coupling gradient. However, if the interstate coupling gradient is zero by definition (for example, for states of different multiplicities), the algorithms can be used with the <tt>coupled_states</tt> option set to false, in which case only the two state gradients are required. The penalty function method never requires the interstate coupling gradient, but is usually less efficient than the other two methods.</p>

<p>Note that the <tt>pf</tt> and <tt>gp</tt> conical intersection optimisation 
algorithms are independent of the choice of <tt>optimiser</tt> in principle 
(in practice <tt>nr</tt> is recommended with <tt>bfgs</tt> updating). 
The Lagrange-Newton method can only be used 
with the Lagrange-Newton optimiser.</p> 



<h3><a name="po" />2.7 Population optimisation options</h3>

<TABLE CELLPADDING=3 BORDER="1">
<TR><TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw1><I>Argument</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw2><I>Argument type</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw3><I>Mandatory</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw4><I>Default</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw5><I>To specify</I></TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">stochastic</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Specifies if stochastic search is to be used.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">genetic</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Specifies if the genetic algorithm is to be used.</TD> 
</TR>
<TR><TD colspan=5><b> General options </b></TD></TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">po_pop_size</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">25</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Working population size.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">po_radius</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">1.0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Maximum displacement when generating initial population.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">po_tolerance_g</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.001</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Threshold for convergence of maximum component of gradient</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">po_maxcycle</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">10000</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Maximum number of optimisation cycles (note <tt>maxcycle</tt> is ignored).</TD> 
</TR>
<TR><TD colspan=5><b> Stochastic search only </b></TD></TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">po_contraction</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.9</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Factor by which search radius is contracted at end of each cycle.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">po_tolerance_r</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">1.0e-8</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Threshold for stopping optimisation based on search radius.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">po_distribution</TD>
<TD ALIGN="LEFT" VALIGN="TOP">keyword</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">force_bias</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Search strategy. Can be <tt>random</tt>, <tt>force_direction_bias</tt> or <tt>force_bias</tt></TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">po_scalefac</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">10.0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Scaling factor for absolute gradient vector in force_bias strategy.</TD> 
</TR>
<TR><TD colspan=5><b> Genetic algorithm only </b></TD></TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">po_init_pop_size</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">2*po_pop_size</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Initial population size from which the lowest energy individuals are taken 
to from the first working population.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">po_reset</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">500</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Number of cycles before resetting population.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">po_mutation_rate</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.15</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Mutation rate (probability of displacement of a coordinate).</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">po_death_rate</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0.5</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Death rate (probability of removing an individual).</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">po_nsave</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">10</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Number of low-energy minima to store.</TD> 
</TR>
</TABLE>

<p>Two population-based search algorithms are implemented in DL-FIND: 
a genetic algorithm and stochastic search. Typically these are used to search 
for global minima on the potential energy surface although the stochastic search 
can also be used with a force bias to find local minima.</p>

<p>The genetic algorithm and stochastic search methods may be run in parallel 
using the <a href="parallel.html#taskfarm">task-farming framework</a>.</p>

<h2><a name="microiter" />3. Microiterative optimisation</h2>

<TABLE CELLPADDING=3 BORDER="1">
<TR><TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw1><I>Argument</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw2><I>Argument type</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw3><I>Mandatory</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw4><I>Default</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw5><I>To specify</I></TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">microiterative</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Specifies if microiterative optimisation is to be used.</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">inner_atoms</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use list_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">none</td> 
<TD ALIGN="LEFT" VALIGN="TOP">List of atoms in the inner region (which should contain all QM atoms).</TD> 
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">inner_residues</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use list_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">none</TD>
<TD ALIGN="LEFT" VALIGN="TOP">List of HDLC residues in the inner region (which should contain all QM atoms). Alternative to <tt>inner_atoms</tt> for HDLC optimisations.
</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">include_res</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Specifies how to handle HDLC residues that cross the boundary between 
inner and outer regions. If false, split any HDLC residues that cross the boundary into two.
If true, include the whole residue in the inner region.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">maxmicrocycle</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">100</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Maximum number of microiterative cycles 
permitted before switching back to macroiterations.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">micro_esp_fit</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">For microiterations, 
approximate the QM region using charges fitted 
to the electrostatic potential. <tt>qm_theory</tt> must be GAMESS-UK, 
TURBOMOLE or MNDO.</TD>
</TR>
</TABLE>

<p>Microiterative optimisation methods improve the efficiency of QM/MM optimisation by 
separating the active atoms into an inner region (which should contain the QM region) 
and an outer environment. After each step taken in the inner region, the environment 
is relaxed completely. The idea is to minimise expensive inner region (macroiterative) 
evaluations at the cost of increasing the number of environment (microiterative) cycles.</p>

<p>Microiterative methods only save overall calculation time if QM region calculations
are not performed during the microiterations. If <tt>micro_esp_fit</tt> is set to true, 
the electrostatic influence of the QM region will be approximated by fitting point charges at the 
QM atom sites to an electrostatic potential generated by the QM code. 
ESP fitting for microiterative optimisation is currently available in the
GAMESS-UK, TURBOMOLE and MNDO interfaces.</p>

<p>Microiterative optimisation has been implemented for minimisation with L-BFGS, transition 
state optimisation with P-RFO and the dimer method, and reaction path optimisation with NEB. 
In all cases the outer environment region is relaxed using L-BFGS. For the transition state 
and reaction path methods, this is equivalent 
to specifying spectator degrees of freedom (setting weights to zero) in a standard optimisation.
This is useful for eliminating complications that can be caused by irrelevant degrees of freedom.
For P-RFO this also means that the Hessian is only calculated over the inner degrees of 
freedom, which can dramatically reduce the cost of the Hessian calculation.</p>

<h2><a name="hessian" />4. Hessian calculations</h2>

<TABLE CELLPADDING=3 BORDER="1">
<TR><TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw1><I>Argument</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw2><I>Argument type</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw3><I>Mandatory</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw4><I>Default</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw5><I>To specify</I></TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">thermal</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">false</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Don't optimise. Calculate the Hessian and 
thermal corrections to the enthalpy and entropy instead.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">temperature</TD>
<TD ALIGN="LEFT" VALIGN="TOP">real</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">300</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Temperature (in Kelvin) at which thermal 
coorections should be calculated.</TD>
</TR>
</TABLE>

<p>DL-FIND may also be used to calculate a finite-difference Hessian and 
thermal corrections to the enthalpy and entropy. In this mode there is no optimisation stage.</p>

<p>If there are no frozen atoms, the rotational and translation modes are projected out 
before determining the vibrational frequencies. If frozen atoms are found this projection is not
applicable but the softest modes can be selectively ignored in the 
thermochemical analysis using the <tt>nzero</tt> option (default: 0).</p>

<p>Note that <tt>thermal</tt> calculations work in mass-weighted coordinates throughout. The 
displacements used and resulting Hessian will therefore not agree exactly with 
the output of a <tt>force</tt> module calculation
(even if the same masses and del value are specified), where the mass-weighted Hessian is 
generated from an initial Cartesian Hessian.</p>

<p>The one-point or two-point finite-difference Hessian may be calculated in parallel
using the <a href="parallel.html#taskfarm">task-farming framework</a>.</p>


<h2><a name="rates" />5. Rate Calculations</h2>

Reaction rates can be calculated in DL-FIND in two different ways: (1)
approximations which only use the energies and Hessians at a minimum
and a saddle point (transition structure) associated to it (i.e. rates
by harmonic transition state theory), and (2)
Instanton rates in which the most-likely tunnelling path is optimised.

<TABLE CELLPADDING=3 BORDER="1">
<TR><TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw1><I>Argument</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw2><I>Argument type</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw3><I>Mandatory</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw4><I>Default</I></TD>
<TD ALIGN="LEFT" VALIGN="TOP" WIDTH=$cw5><I>To specify</I></TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">thermal</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">false</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Don't optimise. Calculate the Hessian and 
thermal corrections to the enthalpy and entropy instead.</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">rate</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">false</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Don't optimise. Read Hessians of the
  reactant state (RS) and the transition state (TS) and calculate
  rates by different methods. No energy or gradient calculations are
  performed.</TD>
</TR>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">qts</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">false</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Search for an instanton path (a quantum transition state)</TD>
</TR>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">qtsrate</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">false</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Calculate the Hessian of an instanton and the
  instanton rate</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">tsplit</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use bool_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">false</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Calculate the tunnelling splitting. Only
  applicable with qtsrate=true</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">inithessian</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">5: skip energy-and-gradient
  calculations, read Hessian from file. 6: Read Hessian of individual
  images from file</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">mass</TD>
<TD ALIGN="LEFT" VALIGN="TOP"><use list_link></TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">masses of the most abundant isotope</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Atomic masses in atomic mass units
  (m(C)=12). One entry per atom. The masses used are printed with list_option=debug</TD>
</TR>
<TR>
<TD ALIGN="LEFT" VALIGN="TOP">nzero</TD>
<TD ALIGN="LEFT" VALIGN="TOP">integer</TD>
<TD ALIGN="LEFT" VALIGN="TOP">no</TD>
<TD ALIGN="LEFT" VALIGN="TOP">0</TD>
<TD ALIGN="LEFT" VALIGN="TOP">Number of modes assumed to have zero vibrational frequency</TD>
</TR>
</TABLE>

<h3><a name="htst" />5.1  Rates Calculated by Harmonic Transition
  State Theory</h3>

 <p> The Hessians at both the reactant state (RS) minimum and the
  transition state (TS) saddle point geometries have to be <a href="#hessian">calculated
  by setting <tt>thermal=true</tt></a>. This writes the files <tt>qts_reactant.txt</tt> and
  <tt>qts_hessian_rs.txt</tt> in case of the reactant, and <tt>qts_ts.txt</tt>
  and <tt>qts_hessian_ts.txt</tt> in case of the transition state. In case of
  a TS, the crossover temperature T<sub>c</sub> is calculated:
  <p> 
  <use equation image=tcrossover.gif 
     latex="$$T_c=\frac{\hbar \omega_b}{2\pi k_B}$$"></p>
  with <use equation image=omegab.gif  latex="$$\omega_b$$"> being the absolute value of the imaginary
   frequency.

  After the Hessians have been calculated, rates are calculated in an
  independent run of dl-find by setting <tt>rate=true</tt>. The files <tt>qts_reactant.txt</tt>,
  <tt>qts_ts.txt</tt>, and <tt>class.in</tt> have to be provided. The latter is an
  input file of the following format:<br> 
  first line: ignored<br>
  second line: number of zero eigenvalues for the reactant and for the TS
  (e.g. <tt>6 6</tt>)<br>
  third line: starting temperature, end temperature, number of
  temperature steps (e.g. <tt>300. 150. 20}</tt>).<br>
  fourth line: <tt>T</tt> if <a href="bimol">bimolecular rates</a> should be calculated.<br>

  The output (stdout or the file <tt>arrhenius</tt>) can directly
  be used in an Arrhenius plot: the first column contains
  1000/T in Kelvin, the next columns contain the
  log<sub>10</sub> of the rates in s<sup>&#8722;1</sup>
  (cm<sup>3</sup>s<sup>&#8722;1</sup> in case
  of <a href="bimol">bimolecular</a> rates) in the following
  approximations:
  <ul>
  <li> Classical rates
  <li> Rates with quantised vibrations (which includes
  the zero-point vibrational energy) 
  <li> Quantised vibrations plus tunnelling correction 
  via the simplified Wigner correction to second order:
  <p> 
  <use equation image=swig2.gif 
     latex="$$\kappa(T)=1+\frac{1}{24} (\beta \hbar \omega_b)^2=1+\frac{1}{24}
    \left(\frac{2\pi T_c}{T}\right)^2, \quad \beta=1/k_BT$$"></p>
   The rate with quantum vibrations is multiplied by the expression above.
  <li> Exact analytical quantum rates for a
  symmetric Eckart barrier fitted to the particular system (barrier
  hight and <use equation image=omegab.gif  latex="$$\omega_b$$">). All degrees of freedom
  perpendicular to the reaction coordinate are approximated as quantum
  mechanical harmonic oscillators.

  <li> For temperatures above the crossover temperature T<sub>c</sub>, the full
  Wigner-corrected rates is also given:
  <p> 
  <use equation image=fwigner.gif 
     latex="$$\kappa(T)=\frac{\beta \hbar \omega_b/2}{\sin(\beta \hbar \omega_b/2)}$$"></p>
  </ul>

  Another file called <tt>arrhenius_polywigner</tt> is written in
  which the rate is calculated by simplified Wigner corrections to second (same
  as above), fourth, sixth, and eighth order. 

  <p>KIEs can be calculated directly by first running DL-FIND with
  <tt>rate=true</tt> using the Hessians for the light
  isotopologue. The file <tt>arrhenius</tt> of this run can be copied to
  the file <tt>rate_H</tt> in the directory where the rate with heavier isotopes is to
  be calculated. There, the same class.in as in the light case (at least the
  same temperature parameters) should be used. The rates obtained with the
  light nuclids are read and the KIE is directly calculated and written to a
  file called <tt>kie</tt>.

  <h3><a name="instanton" />5.2  Rates by Instanton Theory</h3>

  Instanton theory is applicable for temperatures below T<sub>c</sub>
  for each system.

  <p><b>Optimisation of the first instanton starting from the
    classical transition state structure</b> (<tt>qts=true</tt>): The
    file <tt>qts_hessian_ts.txt</tt> has to be renamed
    to <tt>qts_hessian.txt</tt>. All geometrical data are read in
    from <tt>qts_hessian.txt</tt>. However, <tt>coords</tt>
    and <tt>coords2</tt> still have to be provided (for historic
    reasons, number of atoms, ...), but are ignored (as in all
    instanton optimisations and rate calculations). A finite value of
  <tt>distort</tt> specifies how far the images will be spread along the
  unstable mode of the classical TS, see [<a href="#inst_refs">1</a>]. Newton&ndash;Raphson
  optimisation (<tt>optimiser=NR</tt>) is recommended
  [<a href="#inst_refs">1</a>, <a href="#inst_refs">2</a>]. The number of images of the instanton path is
  specified by <tt>nimage</tt>. Instanton search is activated by <tt>qts=true</tt>.
  The NR optimiser was modified to avoid convergence to higher-order saddle
  points [<a href="#inst_refs">2</a>]. This avoids the collapse of the instanton path to the
  classical TS.

  <p>Instanton searches are performed in mass-weighted coordinates with masses
  consistent with atomic units (electron mass, m<sub>e</sub>). That is, the
  mass of a hydrogen atom (<sup>1</sup>H) is 1837.15 m<sub>e</sub>.  This scales all
  distances up by a factor of 42.695 (=(atomic mass
    unit/m<sub>e</sub>)<sup>1/2</sup>) compared to mass-scaled coordinates. Thus, the
  tolerance criterion (<tt>tolerance</tt>) has to be smaller by the same
  factor to achieve equivalent convergence. A tolerance of 10<sup>&#8722;7</sup> (input as
  <tt>1.E-7</tt> in ChemShell) is usually sufficient, but convergence
  even to a tolerance of 10<sup>&#8722;8</sup>  can often still be
  achieved, which results in more reliable rates. Since NR converges quadratically, the more
  stringent tolerance generally does not increase the number of steps
  dramatically.

  <p>If NR (or P-RFO) is used, the updated Hessian will be used to
  calculate a preliminary estimate of the rate (if
  <tt>qts_reactant.txt</tt> is available). In that case,
  <tt>qts_hessian_upd.txt</tt> will be written, which contains only
  the updated Hessian. The file <tt>qts_coords.txt</tt> will be
  written. It acts as input for subsequent recalculation of the
  Hessian and a rate calculation.

  <p>Restarting instanton searches: Proper restart information (check files) is
  not written for the time being. Using NR, a restart is possible, though, by
  renaming <tt>qts_hessian_intermediate.txt</tt> (which is written after
  each step) to <tt>qts_hessian.txt</tt> and starting the simulation
  again. It will start from the Hessian and the geometry after the last full
  set of energies has been obtained.

  <p><b>Instanton rate calculation</b> (<tt>qtsrate=true</tt>): <tt>qts_coords.txt</tt>
  from a previous instanton optimisation is read (<tt>coords</tt> and
  <tt>coords2</tt> are ignored). The temperature is also read from
  <tt>qts_coords.txt</tt>. The rate calculation is selected by
 <tt>qtsrate=true</tt>. Hessians at all images and the
  rate are calculated as described in [<a href="#inst_refs">2</a>]. <tt>qts_hessian.txt</tt> is written, which acts as
  input for subsequent instanton optimisations.

  <p>Restarting of rate calculations is also only possible by using the
  Hessian information written for each image in
  <tt>qts_hessian_imageX.txt</tt>. For these files to be read, set
  <tt>inithessian=6</tt>.

  <p><b>Next instanton optimisation in sequential cooling:</b> Starting from
  <tt>qts_hessian.txt</tt> at a previous (in general higher) temperature,
  another instanton is calculated. <tt>Distort</tt> should be zero, all other
  parameters are the same as above. The number of images may be
  increased. For optimal interpolation, the number of new images
  P<sub>n</sub> should be related to the number of old images P<sub>o</sub>
  by P<sub>n</sub> = k P<sub>o</sub> &#8722;<!-- minus sign --> k + 1
  with k&gt;1 being an integer. This ensures k&#8722;1 new images between each pair
  of old images.

  <p><b>Instanton KIEs</b> can be calculated by starting out from a Hessian
  (<tt>qts_hessian.txt</tt>) for a different isotopologue and
  changing the masses in the input. The Hessian will be re-weighted
  accordingly. The instanton geometry has to be re-optimised. The file
  <tt>qts_reactant.txt</tt> obtained with
  changed masses can not be used. Instead, a file
  <tt>qts_hessian_rs.txt</tt> can be provided. The Hessian of the
  reactant obtained from that file will also be re-weighted.

  <p>In an approximation called fixed-path approximation (FPA) one can keep the instanton geometry fixed
  and just change the masses [<a href="#inst_refs">3</a>]. This is done by calculating an instanton
  rate with <tt>inithessian=5</tt> (read the Hessian from file rather
  than recalculating it) and changing the masses.

<p>For lower temperature (compared to T<sub>c</sub>) the number of
images necessary can be kept at bay by adapting the integration grid
 to the potential energy along the instanton path [<a href="#inst_refs">2</a>]. This only
makes sense if the instanton path has reached the reactant
  minimum. The adaptive grid 
can be achieved by setting <tt>nebk=1</tt> (This is not interpreted as
the NEB force constant, but as a parameter which can vary from zero to
one. One corresponds to a fully adaptive grid).

<p>An instanton rate can be calculated from an existing Hessian (i.e. from the
file <tt>qts_hessian.txt</tt>) without any additional energy
  calculations by setting <tt>inithessian=5</tt>.

<p>Hessians of the individual images can be read in (all or just a part)
from files <tt>qts_hessian_imageX.txt</tt>) by setting <tt>inithessian=6</tt>.

<h3><a name="bimol" />5.3  Bimolecular Rates</h3>

Bimolecular rates are at the moment only implemented for one atom reacting
with a molecule. <tt>qts_reactant.txt</tt> refers to the reactant molecule. I.e. it
has 3 degrees of freedom less than the classical TS. <tt>qts_reactant.txt</tt> has to
be adapted manually: the energy of the incoming atom has to be added to the
third line (which contains the energy of the reactant molecule). Additionally,
the mass on the incoming atom (in atomic mass units) should be appended at the
third line (thus, two real values in the third line).

<p>The relative partition function of the incoming atom will be calculated and
replaces the vibrational partition function for three degrees of freedom. The
rate is internally calculated in atomic units (as is the case for uni-molecular
reactions), but will be converted to molecules cm<sup>3</sup> s<sup>&#8722;1</sup> upon output.

<h3><a name="tsplit" />5.4  Tunnelling Splittings</h3>

Tunnelling splittings of the vibrational ground state level can be
calculated by setting <tt>tsplit=true</tt>. Every time a rate is
calculated, the tunnelling splitting is calculated as well. Tunnelling
splittings only make sense for symmetric molecules and barriers.

<h3><a name="inst_refs" />5.5 References Covering the Implementation
  of Tunnelling Rate Methods</h3>

The following references should be cited if the corresponding methods
are used:
<ul>
<li> [1]  Judith B. Rommel, T. P. M. Goumans, and Johannes
  K&auml;stner, <i>J. Chem. Theory Comput.</i>, 2011, 7,  690. DOI: <a href="http://dx.doi.org/10.1021/ct100658y">10.1021/ct100658y</a>
<li> [2] Judith B. Rommel, Johannes K&auml;stner
  <i>J. Chem. Phys.</i>, 2011, 134, 184107. DOI: <a href="http://dx.doi.org/10.1063/1.3587240">10.1063/1.3587240</a>
<li> [3]  Jan Meisner, Judith B. Rommel and Johannes K&auml;stner
  <i>J. Comput. Chem.</i> 2011, 32, 3456. DOI: <a href="http://onlinelibrary.wiley.com/doi/10.1002/jcc.21930/abstract">10.1002/jcc.21930</a></ul>

<h2><a name="examples" />6. Examples</h2>

<ul>

<li> <b>Minimisation</b> in internal coordinates (DLC):
<pre>
dl-find coords=c \
    theory= mndo: $mndo_args \
    coordinates=dlc \
    result=minimum.c
</pre>

<li> <b>Transition-state search</b> using the P-RFO algorithm (equivalient to "Baker"
in newopt):
<pre>
dl-find coords=c \
    theory= mndo: $mndo_args \
    optimiser=prfo trustradius=const \
    delta=0.01 update_method=bofill maxupdate=50 \
    result=ts_prfo.c result2=tsmode_prfo.c tsrelative=true 
</pre>

<li> Minimisation with internal <b>constraints</b>. Note that internal constraints
(bonds, angles, ...) are only available in DLC and HDLC coordinates. Atoms can
be fixed (frozen) in all coordinate systems using the <tt>active_atoms</tt> or <tt>frozen</tt> options.
<pre>
dl-find coords=c \
    theory= mndo: $mndo_args \
    coordinates=dlc \
    constraints= { { bond 1 2 } { torsion 1 2 3 4 } } 
</pre>

<li><b>Nudged-elastic band</b>: 
It is recommended to use a larger tolerance in NEB calculations than in
normal TS search. The TS can subsequently be refined with the dimer method.
<pre> 
dl-find coords=rs.c coords2= { ts.c ps.c } \
    theory= mndo: $mndo_args \
    coordinates=cartesian \
    optimiser=lbfgs trustradius=const \
    neb=frozen nimage=10 nebk=0.01 \
    tolerance=0.0045 \
    maxcycle=400 maxene=3000 \
    result=ts_neb.c restult2=tsmode_neb.c \
    maxstep=0.9
</pre>

<li> <b>Dimer method</b> for subsequent refinement of the NEB transition state:
<pre>
dl-find coords=ts_neb.c coords2=tsmode_neb.c \
    theory= mndo: $mndo_args \
    coordinates=cartesian \
    optimiser=lbfgs tolerance=0.00045 trustradius=const \
    dimer=true delta=0.01 \
    maxcycle=400 maxene=300 \
    result=ts_dimer.c result2=tsmode_dimer.c \
    maxstep=0.5
</pre>
The transition mode can be visualised by writing a sequence in xyz format:
<pre>
tsmode_xyz coords=ts_dimer.c coords2=tsmode_dimer.c \
    n=20 delta=1. file=tsmode_dimer.xyz
</pre>

<li> To minimise the system starting from a transition state, use the keyword
<tt><b>distort</b></tt>. Note that any coordinate system can be used
(e.g. <tt>coordinates=dlc</tt>) even though the transition mode has been
calculated in cartesians.
<pre>
dl-find coords=ts_dimer.c coords2=tsmode_dimer.c distort=0.1 \
    theory= mndo: $mndo_args \
    coordinates=dlc \
    optimiser=lbfgs trustradius=energy \
    result=minimum.c
</pre>

<li> <b>Conical intersection optimisation</b> using the Lagrange-Newton algorithm:
<pre>
dl-find coords=c \
    theory= mndo: $mndo_args \
    coordinates=cartesian \
    conint=ln state_i=1 state_j=2 \
    optimiser=ln tolerance=0.00045 trustradius=const \
    initial_hessian=diagonal update_method=bfgs maxupdate=1000 \
    maxcycle=400 maxene=300 \
    result=conint.c \
    maxstep=0.1 minstep=0.0
</pre>
<p>Note that <tt>minstep</tt> should be set to zero as the Hessian must always be updated
to ensure that the Lagrange-Newton constraints are imposed. The <tt>ln</tt> optimiser must be 
used with BFGS updates recommended.</p>

<li> A <b>potential-energy-surface scan</b> is described in the <a
href="../tutorial/pe_surface.html#dlf">tutorial</a>. 
</ul>

<file include="tail.htp">
