\documentclass{article}

 \usepackage{textcomp}
 \usepackage{graphicx}
 \usepackage{pifont}
 \usepackage{amsmath}
% \usepackage{floatflt} %floating figures
 \usepackage{varioref}
% non-standard package, use cite instead
% \usepackage{citesort}
\usepackage{cite}
\usepackage{comment}

\usepackage{mathptmx}


\begin{document}

\begin{center}
\LARGE
DL-FIND

\Large
A Geometry Optimiser for Atomistic Simulation Codes

\bigskip
\LARGE
Project Plan and User Manual

\bigskip
\normalsize    
$ $Revision$ $
 
Johannes K\"astner, Thomas Keal, Joanne Carr\\
STFC Daresbury Laboratory,\\
\today
\end{center}

Note: the interface to the calling program is documented in
section~\ref{sec:interface} on page~\pageref{sec:interface}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{General aim and target}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

DL-FIND is a modern and flexible structure optimiser to be included in
electronic structure codes. It provides a stable way to find reaction energy
differences as well as activation barriers. This can be performed starting
from input structures in the region of the reactant and the product.  The
whole process should require as little action from the user as possible.  The
code should parallelise at least the energy and gradient evaluations and
possibly expensive calculations within the optimiser. All units write restart
information at regular intervals to enable a restart of the optimiser.

DL-FIND incorporates an MPI split-communicator taskfarming 
parallelisation framework. This is employed in the parallel optimisers (genetic algorithm 
and stochastic search), finite-difference Hessian evaluation, the 
nudged elastic band method and instanton calculations.

DL-FIND should be cited as: Johannes K\"{a}stner, Joanne M. Carr, Thomas W. Keal, Walter Thiel, Adrian Wander, and Paul Sherwood,
\textit{J.~Phys.~Chem.~A}, 2009, 113 (43), 11856-11865.

A light-weighted description of some capabilities of DL-FIND can be found in
the 2007 issue of Frontiers: 

\texttt{http://www.cse.scitech.ac.uk/about\_us/Frontiers2007/Kaestner - Finding Minima - CSE Frontiers 2007.pdf}

Excited state and parallel optimisation methods are covered in the 2009 Frontiers article: 

\texttt{http://www.cse.scitech.ac.uk/about\_us/Frontiers2009/Keal - Geometry Optimisation - Frontiers2009.pdf}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Functionalities}

\emph{Emphasised:} remains to be implemented

\subsubsection{Coordinate systems:}
\begin{itemize}
\item Cartesians (including frozen atoms and components)
\item Mass-weighted cartesians (including frozen atoms and components)
\item Internals (including all constraints): 
  \begin{itemize}
  \item DLC (delocalised internal coordinates)
  \item DLC-TC (total connection)
  \item HDLC
  \item HDLC-TC
  \end{itemize}
\item \emph{Fractional coordinates and unit cell optimisation?}
\item \emph{Parallel optimisation in internal coordinates?}
\end{itemize}

\subsubsection{Combinations of coordinates (images):}
\begin{itemize}
\item NEB (nudged elastic band)
\item \emph{(Growing) String method}
\item Dimer method \cite{hen99,hey05,kae08}
\item \emph{Replica path method following \cite{woo03}}
\end{itemize}
All of the combinations should work with all versions of coordinate systems.

\subsubsection{Optimisers:}
\begin{itemize}
\item steepest descent
\item conjugate gradient
\item Newton-Raphson/quasi-Newton
\item Damped dynamics
\item L-BFGS
\item P-RFO Hessian update mechanisms: Powell \cite{pow71} and Bofill
  \cite{bof94}. Hessian either by input or by finite-difference. In the latter
  case either in cartesians (then the update also in cartesians, and one can
  output frequencies), or in internals.

  A criterion can be specified with absolute eigenvalues of the hessian below
  that criterion are frozen (considered to be ``soft''). At maximum 6
  eigenmodes are frozen.
\item \emph{Reaction path following / IRC ?}
\end{itemize}

\subsubsection{Microiterative optimisation for QM/MM}
\begin{itemize}
\item Microiterative minimisation with L-BFGS for both regions
\item Microiterative TS search with P-RFO inner region/L-BFGS environment
\item Microiterative TS search with dimer inner region/L-BFGS environment
\item Microiterative NEB with L-BFGS minimised environment
\end{itemize}

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.

Microiterative methods only save overall calculation time if QM region calculations 
are not performed during the microiterations. In the case of electrostatic embedding 
QM/MM calculations with ChemShell, the electrostatic influence of the QM region is 
approximated by fitting point charges at the QM atom sites to an electrostatic 
potential generated by the QM code. 

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.

\subsubsection{Line search algorithms:}
\begin{itemize}
\item Simple scaling of the proposed step (covering the maximum step length)
\item Trust radius based on energy decrease
\item Trust radius based on the projection of the gradient on the step
\item \emph{Trust radius based on the overlap of the lowest eigenmode (for P-RFO)}
\item \emph{full line search ? which algorithm?}
\end{itemize}

\subsubsection{Conical intersection optimisations:}
\begin{itemize}
\item Penalty function
\item Gradient projection
\item Lagrange-Newton
\end{itemize}

\subsubsection{Population-based search:}
\begin{itemize}
\item Random (stochastic) search \cite{brooks57,luusj73}
\item Genetic algorithm \cite{holland75,goldberg89,haupth98}
\end{itemize}

\subsubsection{A global task manager:}

\emph{Should define which methods are used depending on the input. In case a
  method fails, this should be recognised and another method should be tried.}

\subsubsection{Parallelisation:}

The parallel optimisers (stochastic search and genetic algorithm), finite-difference
Hessian evaluation and the nudged elastic band method can be run
in taskfarming mode using MPI (message passing interface), where each taskfarm (workgroup)
calculates the energy and gradient 
for a non-overlapping subset of the total required gradients. The load-balancing 
is static, and the number of workgroups must be a factor of the total number of 
processors for a given job.  If the number of workgroups requested is less than the
total number of processors, then the single-point energy and gradient calculations 
for each individual in a workgroup can also be parallelised.  This is handled by
the program that provides the energy and gradient routines.  Thus, two-level 
parallelisation is possible (as implemented, for example, 
in the task-farming version of ChemShell).  

Wrappers for the required MPI routines are located in the file dlf\_mpi.f90, and the 
corresponding ``dummy'' subroutines are in dlf\_serial.f90.  Which of the corresponding 
.o object files is linked should depend on the build option chosen.  For example, the default build 
of the standalone DL-FIND with its test driver program gives a serial executable.
A \texttt{make parallel} command will produce the parallel-enabled standalone executable, 
Pfind.x.  See makefile and Makefile.standalone for details.

Notes:
\begin{itemize}
\item If DL-FIND is in charge of disentangling the standard output from all the processors, then the 
variable \texttt{keep\_alloutput} in \texttt{dlf\_global\_module.f90} is used.
Rather than being an input option, 
it is hardwired in the code (so that it can be known before any output occurs).
For example, \texttt{keep\_alloutput} is used in \texttt{main.f90} (when a standalone program is made for testing) in the 
interface subroutine \texttt{dlf\_output} as follows: if true, then for processor $n$ ($n/=0$) in the global communicator (\texttt{MPI\_COMM\_WORLD}), 
the file output.proc$<n>$ is opened on unit \texttt{stdout}; 
if false, then for processor $n$ ($n/=0$), /dev/null is ``opened'' on unit \texttt{stdout}.
Output from the rank-zero processor in either case goes to standard out or a named file on unit \texttt{stdout}.

If such I/O issues are dealt with by the main, calling program, then simply use the subroutine \texttt{dlf\_output} to 
change the DL-FIND defaults (\texttt{stdout} and \texttt{stderr} stored in 
\texttt{dlf\_global\_module.f90}) for the required unit numbers, if necessary.
%Default value is .false.

\item Integer variables \texttt{mpi\_rk} and \texttt{mpi\_ik} are declared in \texttt{dlf\_mpi\_module} and set in 
subroutine \texttt{dlf\_mpi\_initialize} to the values of the MPI datatypes 
\texttt{MPI\_Double\_precision} and \texttt{MPI\_Integer}, respectively.
This matches the use of real(rk) in DL-FIND's declarations, where rk = kind(1.d0).  However, problems will arise if 
a compiler flag is used to change the nature of double precision numbers, as the MPI libraries will 
(probably) have been complied with a different specification for double precision.  Could use \texttt{MPI\_Sizeof} and 
\texttt{MPI\_Type\_match\_size} instead if this turns out to be a common problem.

\item Random numbers are an integral part of the parallel optimisers.
Therefore, the seeding of the random-number generator is dealt with in dlf\_parallel\_opt.f90.  Current policy 
is that only the rank-zero processor generates random numbers for the parallel optimisers.  Only the routines that manipulate the population 
need random numbers, so currently only the rank-zero processor does any work in such routines.
\end{itemize}

\subsubsection{Restart mechanism:}

The optimiser should be fully restartable. Status: everything is
restartable, except for the instanton part. There, most information
during the runs is stored as well, but one can not directly start from
the global restart files and continue exactly where the last run has
ended (a practical solution, but not consistent with the rest of DL-FIND)

\subsubsection{Compilation:} The portland compiler version 7.1 and 7.1-1 (and
7.1-2) do not compile the code properly. They do not save the contents of the
hdlc derived type. This is a compiler bug. The portland compilers 6.1 and
7.0-4 do compile it (as well as v9 and later).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Quality Assurance Plan}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

An automated testing system (Buildbot) is used to test DL-FIND as part 
of the ChemShell and GAMESS-UK distributions. 
The ChemShell build of DL-FIND is tested with PGI Fortran, Intel Fortran, g95 and GNU Fortran.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Software Design Plan}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Language: Fortran 95 + TR 15581 (technical report: meaning, allocatable arrays
in derived types are used). Pointers are only to be used where strictly
necessary.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Interface to the Calling Program / API \label{sec:interface}}

DL-FIND is designed as a library to be linked to quantum chemical or MD codes.
However, for testing purposes, a driver module (\texttt{main.f90}) providing some
analytic energy functions is used. DL-FIND is included in GAMESS-UK and
ChemShell.

The routine to be called by the main program is
\texttt{dl\_find(nvarin,nvarin2,nspec,master)}, in the file \texttt{dl\_find.f90}. It
only returns after the complete geometry optimisation.

The main program has to provide the following routines, which are called from
\texttt{dl\_find.f90} and \texttt{dlf\_formstep.f90}:
\begin{itemize}
\item \texttt{dlf\_get\_params(...)} Provide input parameters. The argument
  list is expected to change as DL-FIND is developed. New arguments should
  only be appended to existing ones which makes it possible to keep interfaces
  up to date. The arguments are documented in \texttt{dlf\_global\_module.f90}.
\item \texttt{dlf\_get\_gradient(nvar,coords,energy,gradient,iimage,status)}
  Calculate the energy and gradient at the position coords.
\item \texttt{dlf\_get\_hessian(nvar,coords,hessian,status)} Calculate the
  Hessian at the position coords.
\item \texttt{dlf\_get\_multistate\_gradients(nvar,coords, energy, gradient,
    coupling, needcoupling, iimage, status)} Calculate gradients of a pair of electronic states 
and optionally the interstate coupling between them. Required for conical
    intersection search.
\item \texttt{dlf\_put\_coords(nvar,switch,energy,coords,iam)} Feeds a geometry
  back to the calling program. If switch is 1, coords contains the actual
  geometry. If switch is 2, coords contains the transition mode.
  The presence of iam allows the behaviour on the rank zero processor to be coded 
  differently: for 
  example, turning off the writing of coordinates to a file from all but the
  rank-zero processor.
\item \texttt{dlf\_error()} Error termination. Return to the calling
  program. \texttt{dlf\_error} should not return.
\item \texttt{dlf\_update()} Allows the calling code to update any neighbour
  list, i.e. allows for discontinuities in the potential energy surface. This
  routine is called after a reset of the optimisation algorithm.
\end{itemize}
These may be C routines but have to be Fortran-callable. Examples of interface
routines are available in \texttt{main.f90}, \texttt{dlf.c}, and
\texttt{dlfind\_gamess.m}.

If the taskfarming functionality is required, then the following additions to the 
main program should also be made:
\begin{itemize}
\item call \texttt{dlf\_mpi\_initialize()} to either have DL-FIND set up the global MPI communications, 
      including calling 
      \texttt{MPI\_Init}, or to get the required parameters for the global communications that have 
      already been set up in the 
      main program.  In the former case, this call should be placed before any 
      I/O occurs, as close to the start of the main program as possible, as 
      usual with \texttt{MPI\_Init}.  If the main program sets up the communications, then this 
      call should be made after the rank of each processor and the total number of processors are known 
      (ideally, immediately afterwards).
      This call is an obligatory addition to the main program.
\item call \texttt{dlf\_mpi\_finalize()} to close down MPI before exiting the main program.
      Not necessary if the main program calls \texttt{MPI\_Finalize} already; can be added 
      safely after such a call in the main program (but is redundant).
\item call \texttt{dlf\_output(dlf\_stdout,dlf\_stderr)} passing unit numbers for I/O from 
      the main program to DL-FIND.  The call may be omitted if it is not necessary to change 
      the defaults set in \texttt{dlf\_global\_module.f90}, and if the main program deals with 
      the output from different processors.
\item provide subroutine \texttt{dlf\_output(dlf\_stdout,dlf\_stderr)}, if required, to set 
      the DL-FIND variables stdout 
      and stderr.  If necessary, implement the strategy for dealing with standard output from the 
      different processors here.
\item provide subroutine \texttt{dlf\_put\_procinfo(dlf\_nprocs,dlf\_iam,dlf\_global\_comm)}. 
      Called from \texttt{dlf\_mpi\_initialize} if \texttt{MPI\_Init} was called there.  Passes the 
      total number of processors, 
      the rank of the current process and a variable set to the handle for the 
      global communicator (\texttt{MPI\_COMM\_WORLD)} to the main program. 
\item provide subroutine \texttt{dlf\_get\_procinfo(dlf\_nprocs,dlf\_iam,dlf\_global\_comm)}. 
      Called from \texttt{dlf\_mpi\_initialize} if \texttt{MPI\_Init} had already been called by the 
      main program.  The total number of processors, 
      the rank of the current process and a variable set to the handle for the 
      global communicator (\texttt{MPI\_COMM\_WORLD)} are passed from the main program to 
      DL-FIND. 
\item provide subroutine \texttt{dlf\_put\_taskfarm(dlf\_ntasks,dlf\_nprocs\_per\_task, 
      dlf\_iam\_in\_task,dlf\_mytask,dlf\_task\_comm,dlf\_ax\_tasks\_comm)}.
      Called from \texttt{dlf\_make\_taskfarm} if DL-FIND sets up the split communicators.  
      Passes to the main program the number of taskfarms (workgroups), the number of processors per 
      farm, the rank of 
      the current process in its farm, the rank of the current process's farm, and 
      variables containing the handles for communicators within each farm and for the 
      rank-$n$ processor in each farm.  The main program indicates to DL-FIND which of the 
      two should set up the split communicators via an argument to the general interface routine 
      \texttt{dlf\_get\_params}.  $\texttt{tdlf\_farm}=0$ means the main program does the setup; 
      $\texttt{tdlf\_farm}\neq 0$ means DL-FIND does.
\item provide subroutine \texttt{dlf\_get\_taskfarm(dlf\_ntasks,dlf\_nprocs\_per\_task,
      dlf\_iam\_in\_task,dlf\_mytask,dlf\_task\_comm,dlf\_ax\_tasks\_comm)}.
      Called from dlf\_make\_taskfarm if the main program sets up the split communicators.  
      The number of taskfarms (workgroups), the number of processors per farm, the rank of 
      the current process in its farm, the rank of the current process's farm, and 
      variables containing the handles for communicators within each farm and for the 
      rank-$n$ processor in each farm, are all passed from the main program to DL-FIND.
      The main program indicates to DL-FIND which of the
      two should set up the split communicators via an argument to the general interface routine
      \texttt{dlf\_get\_params}.  $\texttt{tdlf\_farm}=0$ means the main program does the setup;
      $\texttt{tdlf\_farm}\neq 0$ means DL-FIND does.
\item ensure either \texttt{MPI\_Abort} or \texttt{dlf\_mpi\_abort()} is called from the 
      \texttt{dlf\_error()} subroutine.
\end{itemize} 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Program Units and Their Contents:}

\subsubsection{Modules}

The code makes use of four modules that can be used from any subroutine
throughout the DL-FIND: \texttt{dlf\_parameter\_module}, which only provides
the real kind \texttt{rk}. It is located in the file
\texttt{dlf\_stat\_module.f90}.  \texttt{dlf\_checkpoint}: variables and
subroutines for reading and writing checkpoint files. \texttt{dlf\_stat}:
statistics. May be deleted and replaced in the future.  \texttt{dlf\_global}:
contains many global parameters and also arrays (Cartesian coordinates,
internal coordinates, the step, ...). All those parameters are part of a
variable \texttt{glob}, which has a derived type (also defined in the module
\texttt{dlf\_global}). \texttt{pi}, \texttt{stdout}, and \texttt{stderr} are
defined there as well.

All other modules must only be used within the file they are defined! If data
of these module should be provided to other files (units), get- and set
routines must be used!

\subsubsection{Main Units:}

\begin{itemize}
\item Main Unit (\texttt{dl\_find.f90})
\item Convergence tester 
\item Coordinate transform
\item Optimisation algorithms (\texttt{dlf\_formstep.f90} and routines called
  from there)
\item Scalestep: line search and trust radius approaches
\item Utility units
\end{itemize}

\subsubsection{Other details:}

File units that are used longer than for an immediate write ($>1000$):
1001 -- 1001+nimage,max 1050 for xyz of NEB (\texttt{dlf\_neb.f90})
     
When arrays are allocated, the variables \texttt{glob\%storage} and
\texttt{glob\%maxstorage} should be adjusted accordingly to enable control
over the memory usage.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Program Documentation}

The documentation is done using the robodoc tool. Each (important) subroutine
should have an entry marking its connection within the code, the input and
output variables of the global module, and its main purpose. An example file is
\texttt{dlf\_dimer.f90}.

Details of the documentation:

Fortran subroutines are to be documented as functions. The header gives the
unit and the full subroutine name ``NAME'' should not be specified
``SYNOPSIS'' should be the actual code statement beginning the subroutine.

Module header names should be something short and descriptive rather than the
actual Fortran name.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Code fragments obtained from external sources}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Part of \texttt{dlf\_lbfgs.f90} were obtained from
\texttt{http://www.ece.northwestern.edu/\~nocedal/lbfgs.html}.

The HDLC part has been taken from the HDLC optimiser in ChemShell. These are
the files: \texttt{dlf\_hdlc\_constraint.f90} \texttt{dlf\_hdlc\_hdlclib.f90}
\texttt{dlf\_hdlc\_interface.f90} \texttt{dlf\_hdlc\_matrixlib.f90} and
\texttt{dlf\_hdlc\_primitive.f90}. Walter Thiel agreed to make them publically
available (even under GNU license), as long as he gets DL-FIND for the MNDO
code:

\begin{quote}
  [...] I have included the coordinate transformation part of the HDLC
  optimiser into DL-find. I expect HDLCs to work with the dimer method, and
  maybe they even work with NEB. I hope there are no objections from you as
  long as it only goes into ChemShell and gamess.  However, I would like to
  make it available more broadly. It would be interesting to include it into
  other Daresbury codes (I mainly think of DL\_POLY and possibly Crystal -- if
  those people are interested, maybe the solid state people (around Walter
  Temmerman) are also interested).  At a later stage, I think it would also be
  worth making DL-find with the HDLC coordinate transform available under the
  GNU license.  Alternatively, I would have to make a DL-find version that
  runs without HDLCs.
  
  What do think about that? I.e. would you give us (DL) permission to
  distribute part of the HDLCopt code more widely?

  Dear Johannes,
  
  all this is fine with me. It is a good idea to make HDLCopt available to a
  wider audience. In return, I would like to include DL-find in the MNDO
  program (with no restrictions concerning its distribution). I had planned to
  add other optimisers to MNDO anyway, and it would obviously be good not to
  duplicate such work.

  Best wishes and Happy New Year,

  Walter Thiel

\end{quote}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Test Plan}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The driver program will enable tests of all functionalities with analytic
energy function. Should noisy gradients be tested as well, a random part can
be added to the gradient. Test cases from analytic 2-dimensional potentials
showing minima and transition states (e.g. the M\"uller--Brown potential
\cite{mul79}) as well as clusters of Lennard--Jones particles will be used.
The latter can be extended to systems with very many degrees of freedom.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Method documentation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This section documents non-standard implementations in DL-FIND. Published
standard methods are not documented here.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Weights}

A list of weights (one for each atom) can be specified in the input array
coords2. It will be remapped onto weights of each degree of freedom (internal
coordinates) to be optimised.

It can be used to restrict the NEB path to a certain set of atoms, or the
direction of the dimer.

\newpage
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Instanton search and reaction rates}

The code for calculating instantons is not documented in the ChemShell
documentation for the time being. Thus, some user documentation is
summarised here.

\subsubsection{Work flow}

\begin{enumerate}
\item \textbf{Location of a minimum and a saddle point on the potential energy surface}
  (``classical transition state'') associated with the minimum (i.e. no
  barrier between the saddle point and the minimum). Hessians at both of these
  stationary points have to be calculated (ChemShell: \texttt{thermal=true},
  DL-FIND: iopt=11). This writes files \texttt{qts\_reactant.txt} and
  \texttt{qts\_hessian\_rs.txt} in case of the reactant, and \texttt{qts\_ts.txt}
  and \texttt{qts\_hessian\_ts.txt} in case of the transition state. In case of
  a TS, the crossover temperature $T_\mathrm{c}$ is calculated:
  \begin{equation}
    T_\mathrm{c}=\frac{\hbar \omega_b}{2\pi k_\mathrm{B}}
  \end{equation}
   with $\omega_b$ being the absolute value of the imaginary
   frequency.

\item \textbf{Rates without tunnelling:} the files \texttt{qts\_reactant.txt},
  \texttt{qts\_ts.txt}, and \texttt{class.in} have to be provided. The latter is an
  input file of the following format: \\ first line: ignored\\ second line:
  number of zero eigenvalues for the reactant and for the TS
  (e.g. ``\texttt{6~6}'')\\ third line: starting temperature, end temperature,
  number of temperature steps (e.g. \texttt{300.~150.~20}).\\ fourth line:
  ``\texttt{T}'' if bimolecular rates should be calculated, see below
  (\ref{sec:bimol} on page \pageref{sec:bimol})
  
  These files are required to run DL-FIND with \texttt{rate=true} (or
  iopt=13). The output (stdout or the file \texttt{arrhenius}) can directly
  be used in an Arrhenius plot: 1000/$T$ in Kelvin, $\log_{10}$ of the
  classical rates in s$^{-1}$ (cm$^3$s$^{-1}$ in case of bimolecular rates)
  calculated completely classical, with quantised vibrations (which includes
  the zero-point vibrational energy) and including tunnelling approximatively
  via the simplified Wigner correction:
  \begin{equation}
    \kappa(T)=1+\frac{1}{24} (\beta \hbar \omega_b)^2=1+\frac{1}{24}
    \left(\frac{2\pi T_\mathrm{c}}{T}\right)^2 , \quad
    \kappa(T_\mathrm{c})=1+\frac{(2\pi)^2}{24}\approx 2.645
  \end{equation}
  For temperatures above the crossover temperature $T_\mathrm{c}$, the full
  Wigner-corrected rates is also given:
  \begin{equation}
    \kappa(T)=\frac{\beta \hbar \omega_b/2}{\sin(\beta \hbar \omega_b/2)}
  \end{equation}

  In the last column the exact analytical quantum rates for a
  symmetric Eckart barrier fitted to the particular system (barrier
  hight and $\omega_b$) are shown. All degrees of freedom
  perpendicular to the reaction coordinate are approximated as quantum
  mechanical harmonic oscillators.

  KIEs can be calculated directly by first running DL-FIND with
  \texttt{rate=true} (or iopt=13) on the Hessians for the light
  isotopologue. The file \texttt{arrhenius} of this run can be copied to
  \texttt{rate\_H} to 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 nucleids is read and the KIE is directly calculated and written to a
  file called \texttt{kie}.

\item \label{item3} \textbf{Optimisation of the first instanton
    starting from the classical TS:}  The file \texttt{qts\_hessian\_ts.txt}
  has to be renamed to \texttt{qts\_hessian.txt}. All geometrical data
  are read in from \texttt{qts\_hessian.txt}. However, \texttt{coords}
  and \texttt{coords2} 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
  \texttt{distort} specifies how far the images will be spread along
  the unstable mode of the classical TS, see
  \cite{rom11}. Newton--Raphson optimisation (\texttt{optimiser=NR} /
  iopt=20) is recommended \cite{rom11,rom11b}. A QTS search is chosen
  by \texttt{qts=true} (ChemShell) or icoord=190 (DL-FIND).  The NR
  optimiser is modified to avoid convergence to higher-order saddle
  points \cite{rom11b}. This avoids the collapse of the instanton path
  to the classical TS.

  Instanton searches are performed in mass-weighted coordinates with masses
  consistent with atomic units (electron mass, $m_\mathrm{e}$). That is, the
  mass of a hydrogen atom ($^1$H) is 1837.15 $m_\mathrm{e}$.  This scales all
  distances up by a factor of 42.695 ($=(\mbox{atomic mass
    unit}/m_\mathrm{e})^{1/2}$) compared to mass-scaled coordinates. Thus, the
  tolerance criterion (\texttt{tolerance}) has to be smaller by the same
  factor to achieve equivalent convergence. A tolerance of 10$^{-7}$ (input as
  \texttt{1.E-7} in ChemShell) is usually sufficient, a tolerance of 10$^{-8}$
  is also often still possible. Since NR converges quadratically, the more
  stringent tolerance generally does not increase the number of steps
  dramatically.

  If NR (or P-RFO) is used, the updated Hessian will be used to
  calculate a preliminary estimate of the rate (if
  \texttt{qts\_reactant.txt} is available). In that case,
  \texttt{qts\_hessian\_upd.txt} will be written, which contains only
  the updated Hessian. \texttt{qts\_coords.txt} will in any case be
  written. It acts as input for subsequent recalculation of the
  Hessian and a rate calculation.

  Restarting instanton searches: Proper restart information (check files) is
  not written for the time being. Using NR, a restart is possible, though, by
  renaming \texttt{qts\_hessian\_intermediate.txt} (which is written after
  each step) to \texttt{qts\_hessian.txt} and starting the simulation
  again. It will start from the Hessian and the geometry after the last full
  set of energies has been obtained.

\item \textbf{Instanton rate calculation:} \texttt{qts\_coords.txt}
  from a previous instanton optimisation is read (\texttt{coords} and
  \texttt{coords2} are ignored). The temperature is also read from
  \texttt{qts\_coords.txt}. The rate calculation is chosen by
  \texttt{qtsrate=true} or iopt=12. Hessians at all images and the
  rate are calculated as described in
  \cite{rom11b}. \texttt{qts\_hessian.txt} is written, which acts as
  input for subsequent instanton optimisations.

  Restarting of rate calculations is also only possible by using the
  Hessian information written for each image in
  \texttt{qts\_hessian\_imageX.txt}. For these files to be read, set
  \texttt{inithessian=6}.

\item \textbf{Next instanton optimisation in sequential cooling:} Starting from
  \texttt{qts\_hessian.txt} at a previous (in general higher) temperature,
  another instanton is calculated. Distort should be zero, all other
  parameters are the same as in \ref{item3}. The number of images may be
  increased. For optimal interpolation, the number of new images
  $P_\mathrm{n}$ should be related to the number of old images $P_\mathrm{o}$
  by:
  \begin{equation}
    P_\mathrm{n}=k\ P_\mathrm{o} - k + 1
  \end{equation}
  with $k>1$ being an integer. This ensures $k-1$ new images between each pair
  of old images.

\item \textbf{Instanton KIEs} can be calculated by starting out from a Hessian
  (\texttt{qts\_hessian.txt}) 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
  \texttt{qts\_reactant.txt} obtained with
  changed masses can not be used. Instead, a file
  \texttt{qts\_hessian\_rs.txt} (which includes the masses, so from a
  recent version of DL-FIND) can be provided. The Hessian of the
  reactant obtained from that file will also be re-weighted.

  In an approximation (FPA) one can keep the instanton geometry fixed
  and just change the masses \cite{mei11}. This is done by calculating an instanton
  rate with \texttt{inithessian=5} (read the Hessian from file rather
  than recalculating it) and changing the masses.
\end{enumerate}

For lower temperature (compared to $T_\mathrm{c}$) the number of
images necessary can be kept at bay by adapting the integration grid
(dtau) to the potential energy along the instanton path \cite{rom11b}. This only
makes sense if the instanton path has reached the reactant minimum. It
can be achieved by setting \texttt{nebk=1} (this is not interpreted as
the NEB force constant, but as a parameter which can vary from 0 to
1. One corresponds to a fully adaptive grid).

If KIEs should be calculated with \texttt{rate=true} (or iopt=13), it is
not necessary to recalculate the Hessian for the heavier isotopologue. The
Hessians for the reactant and product are read in. If the masses provided via
the calling code (ChemShell or main.f90) are different from the ones in the
Hessian files, the Hessians with the new masses will be calculated.

In older versions of DL-FIND the masses were not written into the file
\texttt{qts\_hessian.txt}. The mass is needed there for calculating KIEs
(unless one wants to recalculate the whole Hessian), however. A workaround is:
delete everything below the second line in \texttt{qts\_reactant.txt}. Rename
the Hessian file for which masses are needed to \texttt{qts\_hessian\_rs.txt},
and run DL-FIND with \texttt{rate=true} (or iopt=13). It is important that the
masses provided to DL-FIND by the calling code are the same as the ones used
to calculate the Hessian. Then, DL-FIND will write a file
\texttt{qts\_hessian\_rs\_mass.txt} with the masses in, which can be used as
\texttt{qts\_hessian.txt} or \texttt{qts\_hessian\_rs.txt} in subsequent
calculations.

An instanton rate can be calculated from an existing Hessian (i.e. from the
file \texttt{qts\_hessian.txt}) by setting \texttt{inithessian=5}.

Hessians of the individual images can be read in (all or just a part)
from files \texttt{qts\_hessian\_imageX.txt}) by setting \texttt{inithessian=6}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Bimolecular rates \label{sec:bimol}}

Bimolecular rates are at the moment implemented in two ways. The simple one is
for one atom reacting with a molecule. \texttt{qts\_reactant.txt} refers to
the reactant molecule. I.e. it has 3 degrees of freedom less than the
classical TS. \texttt{qts\_reactant.txt} 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).

The relative translational 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$^3$
s$^{-1}$ upon output.

For two molecules reacting with each other with more than one atom in each, a
second value has to be added to the third line of \texttt{qts\_reactant.txt}
as well. However, any negative value will do. It is only used as a label. The
Hessians of the reactants are read from the files
\texttt{qts\_hessian\_rs.txt} and \texttt{qts\_hessian\_rs2.txt} which are
obtained from previous DL-FIND runs. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Tunnelling splittings}

Tunnelling splittings of the vibrational ground state level following
\cite{ric11} can be calculated by setting \texttt{tsplit=true} in ChemShell or
qtsflag=1 in DL-FIND. Every time a rate is calculated, the tunnelling splitting
is calculated as well. Tunnelling splittings only make sense for
symmetric molecules and barriers.

\newpage
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Documentation of the Input Options -- User Documentation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

System documentation is available in the source code after the subroutine
headers. The input options will be explained here (User Documentation).

\paragraph{Print level (Variable \texttt{printl}):}
\begin{description}
\item[0] no printout
\item[2] print something
\item[4] be verbose
\item[6] debug
\end{description}

\paragraph{Type of coordinate system (Variable \texttt{icoord}):} 
\begin{description}
\item ``unit place'' means \texttt{icoord} modulo 10
\item[0--9] The whole system is to be treated as one image 
\item[Unit place 0] Cartesians
\item[Unit place 1] HDLC - internals \cite{bil00a}
\item[Unit place 2] HDLC - TC \cite{bil00a}
\item[Unit place 3] DLC - internals \cite{bil00a}
\item[Unit place 4] DLC - TC \cite{bil00a}
\item[Unit place 5] Mass-weighted Cartesians ($\sqrt{m}x$) -- will be deleted:
  use glob\%massweight for that now!
\item[1X] Lagrange-Newton conical intersection search, with two extra 
coordinates corresponding to the gradient difference vector and 
interstate coupling gradient constraints.
\item[10X]  NEB with endpoints free
\item[11X] NEB with endpoints moving only perpendicular to their tangent
  direction
\item[12X] NEB with frozen endpoints.
\item[13X]  NEB with endpoints free. Only initialisation in coordiates X,
  optimisation in cartesians.
\item[14X] NEB with endpoints moving only perpendicular to their tangent
  direction. Only initialisation in coordiates X,
  optimisation in cartesians.
\item[15X] NEB with frozen endpoints. Only initialisation in coordiates X,
  optimisation in cartesians.
  
  The NEB version implemented is the ``improved-tangent'' NEB \cite{hen00a}
  (also called ``upwind scheme'') with a climbing image.

\item[190] Quantum transiton state search
  
\item[20X] Dimer method \cite{hen99,hey05}. Translation and rotation of the
  dimer are covered by the optimiser specified trough \texttt{iopt}. Requires
  two energy evaluation per iteration.
\item[21X] Dimer method. Rotation of the dimer is done by a line search within
  the dimer module, two energy calculations are used per rotation. Requires at
  least two energy evaluation per iteration.
\item[22X] Dimer method. Rotation of the dimer is done by a line search within
  the dimer module, one energy calculation is done per iteration, the other
  one is interpolated. Requires at least two energy evaluation per iteration.
  
\item[30X] ``Chain'' search (similar to NEB, but path is expanded in arbitrary basis functions). Better name required.

  In all dimer versions: If a second set of coordinates is provided, it
  determines the dimer direction, if not, the dimer direction is randomised.
\end{description}

\paragraph{Multistate calculations (Variable \texttt{imultistate}):}
\begin{description}
\item[0] Single state calculation (default)
\item[1] Conical intersection optimisation (penalty function algorithm).
\item[2] Conical intersection optimisation (gradient projection algorithm).
\item[3] Conical intersection optimisation (Lagrange-Newton algorithm).
\end{description}

\paragraph{Type of optimisation algorithm (Variable \texttt{iopt}):}
\begin{description}
\item[0] Steepest descent
\item[1] Conjugate gradient following Polak--Ribi\`ere \cite{pol69} (with automatic
  restarts based on the criterion by Powell and Beale) that is not coded
  properly ...
\item[2] Conjugate gradient following Polak--Ribi\`ere \cite{pol69} with CG
  restart every 10 steps (hardcoded at the moment)
\item[3] L-BFGS \cite{liu89,noc80}
\item[9] Test delta for finite-difference in gradients (19 energy and gradient evaluations)
\item[10] P-RFO \cite{cer81,sim83,ban85,bak86} A switching mechanism for the
  mode to be followed is included, but does not seem to help in any of the
  cases I tried sofar.
\item[11] Just calculate the Hessian and do a thermal analysis (harmonic
  approximation for entropy, ...)
\item[12] Calculate the Hessians of all images and the qTS rate (if inithessian=5: just the rate, read the Hessians)
\item[13] Rate without tunneling (only with Wigner correction)
\item[20] Newton--Raphson/quasi-Newton
\item[30] Damped dynamics using the variables timestep, fric0, fricfac, and
  fricp. The frictions are defined that 0 corresponds to free (undamped)
  dynamics, and 1 corresponds to steepest descent.
\item[51] Random (stochastic) search \cite{brooks57,luusj73}, using the variables po\_pop\_size, po\_radius,
            po\_contraction, po\_tolerance\_r, po\_tolerance\_g, 
            po\_distribution, po\_maxcycle, po\_scalefac.
\item[52] Genetic algorithm \cite{holland75,goldberg89,haupth98}, using the variables po\_pop\_size, po\_radius,
            po\_tolerance\_g, po\_maxcycle, po\_init\_pop\_size, po\_reset, 
            po\_mutation\_rate, po\_death\_rate, po\_nsave
\end{description}

\paragraph{Type of line search or trust radius (Variable \texttt{iline}):}
\begin{description}
\item[0]   simple scaling of the proposed step, taking maxstep into account
\item[1]   Trust radius based on energy as acceptance criterion (recommended
  for L-BFGS optimisation)
\item[2]   Trust radius based on gradient as acceptance criterion (recommended
  for CG optimisation)
\item[3]   Hard-core line search. Does not work at the moment...
\end{description}

\paragraph{Type of initial Hessian (Variable \texttt{inithessian}):}
\begin{description}
\item[0]   Calculate externally using \texttt{dlf\_get\_hessian}. Defaults to two point finite
  difference if an external Hessian is unavailable.
\item[1]   Build by one point finite difference of the gradient
\item[2]   Build by two point finite difference of the gradient
\item[3]   Build a diagonal Hessian with a single one point finite difference
\item[4]   Set the Hessian to be an identity matrix
\item[5]   Only for instanton calculations: read the Hessian from \texttt{qts\_hessian.txt}
\item[6]   Only for instanton calculations: read the Hessian from files
  for each image\texttt{qts\_hessian\_imageX.txt}
\end{description}

\paragraph{Hessian update mechanism (Variable \texttt{update}):}
\begin{description}
\item[0]   No update. Always recalculate the Hessian
\item[1]   Powell update \cite{pow71}
\item[2]   Bofill update \cite{bof94}
\item[3]   BFGS update
\end{description}

\paragraph{Fragment and frozen atom information (Variable \texttt{spec}):}
The array \texttt{spec} has an entry for each atom. Meaning:
\begin{description}
\item[$>$0] Fragment (residue) number the atom belongs to. Active in the
  optimisation.
\item[0]   Active, and treated in Cartesian coordinates, even if other atoms
  in the system are covered by (H)DLCs.
\item[-1]  Frozen 
\item[-2]  x-component frozen (Cartesians only)
\item[-3]  y-component frozen (Cartesians only)
\item[-4]  z-component frozen (Cartesians only)
\item[-23]  x and y-components frozen (Cartesians only)
\item[-24]  x and z-components frozen (Cartesians only)
\item[-34]  y and z-components frozen (Cartesians only)
\end{description}
Atoms with \texttt{spec} $<-1$ will be completely frozen if used in HDLCs. One
may use Cartesian constraints there to specify frozen components. Atoms with
\texttt{spec} $>0$ will be free in Cartesian coordinats.

After this array, \texttt{spec} also contains the following information:
\begin{description}
\item[nz]      entries of nuclear charges (same order as coords)
\item[5*ncons] entries of constraints (typ, atom1,atom2, atom3, atom4)
\item[2*nconn] entries of connections (atom1 atom2)
\item[nat]     entries of microiterative region specification
\end{description}
i.e. \texttt{nspec= nat + nz + 5*ncons + 2*nconn + nat}

\paragraph{Microiterative optimisation (Variable \texttt{imicroiter}):}
\begin{description}
\item[0]   Standard (non-microiterative) optimisation
\item[1]   Microiterative optimisation
\end{description}
Note imicroiter is also used internally to keep track of whether the optimisation 
is in a macroiterative (imicroiter=1) or microiterative (imcroiter=2) loop.

The inner region specification for microiterative optimisation is 
part of the \texttt{spec} array with a section of length \texttt{nat}. A 0 entry signifies 
a standard optimisation or outer region as appropriate, while 1 signifies 
an inner region atom.

\paragraph{Other parameters}
There are numerous other paramers that can be set via the routine
\texttt{dlf\_get\_params}. These are at the moment explained in the global
module, \texttt{dlf\_global\_module.f90}.

\subsubsection{Restarting}

\begin{description}
\item[\texttt{dump}] after how many energy and gradient evaluations is a
  restart file (checkpoint file) to be written? Default: 0 (never).
\item[\texttt{restart}] 0: new run (default), 1: read all checkpoint files and start
  from those.
\end{description}

It is only possible to restart a job with most parameters equal to those in
the checkpoint file. Exeptions (input parameters that may be different from
the ones in the previous run, i.e. that are not written to the checkpoint
file): \texttt{maxcycle}. If different parameters are to be used, the latest
geometry should be used and the optimiser should be started from scratch
(\texttt{restart=0}).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Collection ...}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Force weighted intenal coordinates: have a look at
http://www.molpro.net/molpro-user/archive/all/msg00071.html

At the moment, I am experiencing problems with NEB in internal coordinates.
While the structure may change continously, the internals may not. Consider
H$_2$CO separation: it has two dihedrals. In the bound configuration, both are
180\textdegree{}. In the dissociated configuration, one is 0 and one is 180,
as H$_2$ moves to one side of CO. How to deal with that? 

One possible solution: Impropers are automatically put onto atoms with nearly
planar configuration (what if more than 3 connections?). If an improper is
placed on an atom, this one should not be in the middle of a dihedral that
include these atoms. This, however, had be removed again as it leads to an
underdetermined systems when not all atoms next to the planar atom are
monovalent.

svn properties:

\verb'svn propset svn:keywords "URL Author Date Rev Id" dlf_util.f90'

%\bibliographystyle{myunsrt}
%\addcontentsline{toc}{chapter}{Bibliography}
%\bibliography{text}

\begin{thebibliography}{10}

\bibitem{hen99}
G.~Henkelman and H.~J{\'o}nsson:
\newblock A dimer method for finding saddle points on high dimensional
  potential surfaces using only first derivatives.
\newblock {\em J. Chem. Phys.\/}~{\bf 111}, 7010  (1999).

\bibitem{hey05}
A.~Heyden, A.~T. Bell, and F.~J. Keil:
\newblock Efficient methods for finding transition states in chemical
  reactions: Comparison of improved dimer method and partitioned rational
  function optimization method.
\newblock {\em J. Chem. Phys.\/}~{\bf 123}, 224101  (2005).

\bibitem{kae08}
J.~K\"astner and P.~Sherwood:
\newblock Superlinearly converging dimer method for transition state search.
\newblock {\em J. Chem. Phys.\/}~{\bf 128}, 014106  (2008).

\bibitem{woo03}
H.~L. Woodcock, M.~Hodo{\v s}{\v c}ek, P.~Sherwood, Y.~S. Lee, H.~F.
  Schaefer~III, and B.~R. Brooks:
\newblock Exploring the quantum mechanical/molecular mechanical replica path
  method: a pathway optimization of the chorismate to prephenate Claisen
  rearrangement catalyzed by chorismate mutase.
\newblock {\em Theor. Chem. Acc.\/}~{\bf 109}, 140  (2003).

\bibitem{pow71}
M.~J.~D. Powell: .
\newblock {\em Math. Prog.\/}~{\bf 26}, 1  (1971).

\bibitem{bof94}
J.~M. Bofill:
\newblock Updated Hessian matrix and the restricted step method for locating
  transition structures.
\newblock {\em J. Comput. Chem.\/}~{\bf 15}, 1  (1994).

\bibitem{brooks57}
S.~H.~Brooks:
\newblock A discussion of random methods for seeking maxima.
\newblock {\em Operations Research}~\textbf{6}, 244 (1957).

\bibitem{luusj73}
R.~Luus, and T.~H.~I.~Jaakola:
\newblock Optimization by Direct Search and Systematic Reduction of the Size of Search Region.
\newblock {\em AIChE Journal}~\textbf{19}, (1973).

\bibitem{holland75}
J.~H.~Holland:
\newblock Adaptation in Natural and Artificial Systems.
\newblock University of Michigan Press, Ann
  Arbor, (1975).

\bibitem{goldberg89}
D.~Goldberg:
\newblock Genetic Algorithms in Search, Optimization, and Machine Learning.
\newblock Addison-Wesley, New York, (1989).

\bibitem{haupth98}
R.~L.~Haupt, and S.~E.~Haupt:
\newblock Practical genetic algorithms.
\newblock Wiley, New York, (1998).

\bibitem{mul79}
K.~M{\"u}ller and L.~D. Brown:
\newblock Location of saddle points and minimum energy paths by a constrained
  simplex optimization procedure.
\newblock {\em Theor. Chim. Acta\/}~{\bf 53}, 75  (1979).

\bibitem{rom11}
J.~B. Rommel, T.~P.~M. Goumans, and J.~K\"astner:
\newblock Locating instantons in many degrees of freedom.
\newblock {\em J. Chem. Theory Comput.\/}~{\bf 7}, 690  (2011).

\bibitem{rom11b}
J.~B. Rommel, J.~K\"astner:
\newblock Adaptive Integration Grids in Instanton Theory Improve the Numerical Accuracy at Low Temperature.
\newblock {\em  J. Chem. Phys.\/}~{\bf 134}, 184107 (2011).

\bibitem{mei11}
J.~Meisner, J.~B. Rommel, J.~K\"astner:
\newblock Kinetic Isotope Effects Calculated with the Instanton Method 
\newblock {\em  J. Comput. Chem.\/}~{\bf 32}, 3456 (2011).

\bibitem{ric11}
J.~O. Richardson and S.~C. Althorpe:
\newblock Ring-polymer instanton method for calculating tunneling splittings.
\newblock {\em J. Chem. Phys.\/}~{\bf 134}, 054109  (2011).

\bibitem{bil00a}
S.~R. Billeter, A.~J. Turner, and W.~Thiel:
\newblock Linear scaling geometry optimization and transition state search in
  hybrid delocalised internal coordinates.
\newblock {\em Phys. Chem. Chem. Phys.\/}~{\bf 2}, 2177  (2000).

\bibitem{hen00a}
G.~Henkelman and H.~J{\'o}nsson:
\newblock Improved tangent estimate in the nudged elastic band method for
  finding minimum energy paths and saddle points.
\newblock {\em J. Chem. Phys.\/}~{\bf 113}, 9978  (2000).

\bibitem{pol69}
E.~Polak and G.~Ribi\`ere:
\newblock Note sur la convergence de directions conjug\'ees.
\newblock {\em Ref. Fra. Inf. Rech. Op.\/}~{\bf 3}, 35  (1969).

\bibitem{liu89}
D.~C. Liu and J.~Nocedal:
\newblock On the limited memory BFGS method for large scale optimization.
\newblock {\em Math. Program.\/}~{\bf 45}, 503  (1989).

\bibitem{noc80}
J.~Nocedal:
\newblock Updating quasi-Newton matrices with limited storage.
\newblock {\em Math. Comp.\/}~{\bf 35}, 773  (1980).

\bibitem{cer81}
C.~J. Cerjan and W.~H. Miller:
\newblock On finding transition states.
\newblock {\em J. Chem. Phys.\/}~{\bf 75}, 2800  (1981).

\bibitem{sim83}
J.~Simons, P.~J{\o}rgensen, H.~Taylor, and J.~Ozment:
\newblock Walking on potential energy surfaces.
\newblock {\em J. Phys. Chem.\/}~{\bf 87}, 2745  (1983).

\bibitem{ban85}
A.~Banerjee, N.~Adams, J.~Simons, and R.~Shepard:
\newblock Search for stationary points on surfaces.
\newblock {\em J. Phys. Chem.\/}~{\bf 89}, 52  (1985).

\bibitem{bak86}
J.~Baker:
\newblock An algorithm for the location of transition states.
\newblock {\em J. Comput. Chem.\/}~{\bf 7}, 385  (1986).

\end{thebibliography}

\end{document}


