Geometry optimization
*********************

A geometry optimization is the process of changing the system's
geometry (the nuclear coordinates and potentially the lattice vectors)
to minimize the total energy of the system. This is typically a local
optimization, i.e. the optimization converges to the next local
minimum on the potential energy surface (PES), given the initial
system geometry specified in the "System" block. In other words: The
geometry optimizer moves "downhill" on the PES into the local minimum.

See also:

  Examples and silicon lattice optimization and phonons tutorial

Geometry optimizations are performed by selecting them as the "Task".
The details of the optimization can be configured in the corresponding
block:

   Task GeometryOptimization

   GeometryOptimization
      Convergence
         Energy float
         Gradients float
         Quality [VeryBasic | Basic | Normal | Good | VeryGood | Custom]
         Step float
         StressEnergyPerAtom float
      End
      MaxIterations integer
      CalcPropertiesOnlyIfConverged Yes/No
      OptimizeLattice Yes/No
      KeepIntermediateResults Yes/No
      PretendConverged Yes/No
      MaxRestarts integer
      RestartDisplacement float
   End


Convergence criteria
====================

"GeometryOptimization"
   Type:
      Block

   Description:
      Configures details of the geometry optimization and transition
      state searches.

   "Convergence"
      Type:
         Block

      Description:
         Convergence is monitored for up to 4 quantities: the energy
         change, the Cartesian gradients, the Cartesian step size, and
         for lattice optimizations the stress energy per atom.
         Convergence criteria can be specified separately for each of
         these items.

      "Energy"
         Type:
            Float

         Default value:
            1e-05

         Unit:
            Hartree

         Value Range:
            value > 0

         GUI name:
            Energy convergence

         Description:
            The criterion for changes in the energy. The energy is
            considered converged when the change in energy is smaller
            than this threshold times the number of atoms.

      "Gradients"
         Type:
            Float

         Default value:
            0.001

         Unit:
            Hartree/Angstrom

         Value Range:
            value > 0

         GUI name:
            Gradient convergence

         Description:
            Threshold for nuclear gradients.

      "Quality"
         Type:
            Multiple Choice

         Default value:
            Custom

         Options:
            [VeryBasic, Basic, Normal, Good, VeryGood, Custom]

         GUI name:
            Convergence

         Description:
            A quick way to change convergence thresholds: 'Good' will
            reduce all thresholds by an order of magnitude from their
            default value. 'VeryGood' will tighten them by two orders
            of magnitude. 'Basic' and 'VeryBasic' will increase the
            thresholds by one or two orders of magnitude respectively.

      "Step"
         Type:
            Float

         Default value:
            0.01

         Unit:
            Angstrom

         Value Range:
            value > 0

         GUI name:
            Step convergence

         Description:
            The maximum Cartesian step allowed for a converged
            geometry.

      "StressEnergyPerAtom"
         Type:
            Float

         Default value:
            0.0005

         Unit:
            Hartree

         Value Range:
            value > 0

         Description:
            Threshold used when optimizing the lattice vectors. The
            stress is considered 'converged' when the maximum value of
            stress_tensor * cell_volume / number_of_atoms is smaller
            than this threshold (for 2D and 1D systems, the
            cell_volume is replaced by the cell_area and cell_length
            respectively).

A geometry optimization is considered converged when all the following
criteria are met:

1. The difference between the bond energy at the current geometry and
   at the previous geometry step is smaller than "Convergence%Energy"
   times the number of atoms in the system.

2. The maximum Cartesian nuclear gradient is smaller than
   "Convergence%Gradient".

3. The root mean square (RMS) of the Cartesian nuclear gradients is
   smaller than 2/3 "Convergence%Gradient".

4. The maximum Cartesian step is smaller than "Convergence%Step".

5. The root mean square (RMS) of the Cartesian steps is smaller than
   2/3 "Convergence%Step".

**Note**: If the maximum and RMS gradients are 10 times smaller than
the convergence criterion, then criteria 4 and 5 are ignored.

The "Convergence%Quality" sets the criteria as follows:

+----------------------+----------------------+----------------------+----------------------+----------------------+
| Quality              | Energy (Ha)          | Gradients (Ha/Å)     | Step (Å)             | StressEnergyPerAtom  |
|                      |                      |                      |                      | (Ha)                 |
|======================|======================|======================|======================|======================|
| VeryBasic            | 10⁻³                 | 10⁻¹                 | 1                    | 5×10⁻²               |
+----------------------+----------------------+----------------------+----------------------+----------------------+
| Basic                | 10⁻⁴                 | 10⁻²                 | 0.1                  | 5×10⁻³               |
+----------------------+----------------------+----------------------+----------------------+----------------------+
| Normal               | 10⁻⁵                 | 10⁻³                 | 0.01                 | 5×10⁻⁴               |
+----------------------+----------------------+----------------------+----------------------+----------------------+
| Good                 | 10⁻⁶                 | 10⁻⁴                 | 0.001                | 5×10⁻⁵               |
+----------------------+----------------------+----------------------+----------------------+----------------------+
| VeryGood             | 10⁻⁷                 | 10⁻⁵                 | 0.0001               | 5×10⁻⁶               |
+----------------------+----------------------+----------------------+----------------------+----------------------+

Some remarks on the choice of the convergence thresholds:

* Molecules may differ very much in the stiffness around the energy
  minimum. Using the standard convergence thresholds without second
  thought is therefore not recommended. Strict criteria may require a
  large number of steps, while a loose threshold may yield geometries
  that are far from the minimum (with respect to atom-atom distances,
  bond-angles etc...) even when the total energy of the molecule might
  be very close to the value at the minimum. It is good practice to
  consider first what the objectives of the calculation are. The
  default settings in AMS are intended to be reasonable for most
  applications, but inevitably situations may arise where they are
  inadequate.

* The convergence threshold for the coordinates ("Convergence%Step")
  is not a reliable measure for the precision of the final
  coordinates. Usually it yields a reasonable estimate (order of
  magnitude), but to get accurate results one should tighten the
  criterion on the gradients, rather than on the steps (coordinates).
  (The reason for this is that with the Quasi-Newton based optimizers
  the estimated uncertainty in the coordinates is related to the used
  Hessian, which is updated during the optimization. Quite often it
  stays rather far from an accurate representation of the true
  Hessian. This does usually not prevent the program from converging
  nicely, but it does imply a possibly incorrect calculation of the
  uncertainty in the coordinates.)

* Note that tight convergence criteria for the geometry optimization
  require accurate and noise-free gradients from the engine. For some
  engines this might mean that their numerical accuracy has to be
  increased for geometry optimization with tight convergence criteria,
  see e.g. the "NumericalQuality" keyword in the BAND manual.

The maximum number of geometry iterations allowed to locate the
desired structure is specified with the "MaxIterations" keyword:

"GeometryOptimization"
   "MaxIterations"
      Type:
         Integer

      Value Range:
         value >= 0

      GUI name:
         Maximum number of iterations

      Description:
         The maximum number of geometry iterations allowed to converge
         to the desired structure.

   "CalcPropertiesOnlyIfConverged"
      Type:
         Bool

      Default value:
         Yes

      Description:
         Compute the properties requested in the 'Properties' block,
         e.g. Frequencies or Phonons, only if the optimization (or
         transition state search) converged. If False, the properties
         will be computed even if the optimization did not converge.

   "PretendConverged"
      Type:
         Bool

      Default value:
         No

      Description:
         Normally a non-converged geometry optimization is considered
         an error. If this keyword is set to True, the optimizer will
         only produce a warning and still claim that the optimization
         is converged. (This is mostly useful for scripting
         applications, where one might want to consider non-converged
         optimizations still successful jobs.)

If the geometry optimization does not converge within this many steps
it is considered failed and the iteration aborted, i.e. PES point
properties block will not be calculated at the last geometry. The
default maximum number of steps is chosen automatically based on the
used optimizer and the number of degrees of freedom to be optimized.
The default is a fairly large number already, so if the geometry has
not converged (at least to a reasonable extent) within that many
iterations you should step back and consider the underlying cause
rather than simply increase the allowed number of iterations and try
again.


Automatic restarts
==================

While a geometry optimization aims to find a (local) PES minimum, it
may occur that it ends up finding a saddle point instead. The
"PESPointCharacter" property keyword can be used to quickly calculate
the lowest few Hessian eigenvalues to determine what kind of
stationary PES point the optimization found. More information on this
feature can be found on its Documentation Page. Since AMS2022.1,
geometry optimizations with enabled PES point characterization can
automatically restart when a transition state (or higher order saddle
point) is found: the geometry is distorted along the lowest frequency
mode and the optimizer run again. The applied distortion is often
symmetry breaking in this case, so this automatic restarting is only
enabled if the system does not have any symmetry operators or the use
of symmetry is explicitly disabled using the "UseSymmetry" keyword.
Furthermore the automatic restarting must be explicitly enabled by
setting the "MaxRestarts" option to a value >0. Of course the PES
point characterization needs to be enabled too:

   GeometryOptimization
      MaxRestarts 5
   End

   UseSymmetry False

   Properties
      PESPointCharacter True
   End

Details of the automatic restarting can be configured with the
following keywords:

"GeometryOptimization"
   "MaxRestarts"
      Type:
         Integer

      Default value:
         0

      Description:
         If a geometry optimization of a system with no symmetry
         operators (or with explicitly disabled symmetry: "UseSymmetry
         False") and enabled PES point characterization converges to a
         transition state (or higher order saddle point), it can be
         restarted automatically after a small displacement along the
         imaginary vibrational mode. In case the restarted
         optimization again does not find a minimum, this can happen
         multiple times in succession. This keyword sets the maximum
         number of restarts. The default value is 0, so the automatic
         restarting is disabled by default.

   "RestartDisplacement"
      Type:
         Float

      Default value:
         0.05

      Unit:
         Angstrom

      Description:
         If a geometry optimization of a system with no symmetry
         operators (or with explicitly disabled symmetry: "UseSymmetry
         False") and enabled PES point characterization converges to a
         transition state (or higher order saddle point), it can be
         restarted automatically after a small displacement along the
         imaginary vibrational mode. This keywords sets the size of
         the displacement for the furthest moving atom.


Lattice optimization
====================

For periodic systems the lattice degrees of freedom can be optimized
in addition to the nuclear positions.

"GeometryOptimization"
   "OptimizeLattice"
      Type:
         Bool

      Default value:
         No

      Description:
         Whether to also optimize the lattice for periodic structures.
         This is currently supported with the Quasi-Newton, FIRE, and
         L-BFGS optimizers.

See also Constrained optimization for constrained lattice
optimizations.


Keep all results files
======================

The "GeometryOptimization" block also contains some technical options:

"GeometryOptimization"
   "KeepIntermediateResults"
      Type:
         Bool

      Default value:
         No

      Description:
         Whether the full engine result files of all intermediate
         steps are stored on disk. By default only the last step is
         kept, and only if the geometry optimization converged. This
         can easily lead to huge amounts of data being stored on disk,
         but it can sometimes be convenient to closely monitor a
         tricky optimization, e.g. excited state optimizations going
         through conical intersections, etc. ...


Constrained optimization
========================

The AMS driver also allows constrained optimizations, where a number
of specified degrees of freedom are fixed to particular values.

See also: Example demonstrating all supported constraints

The desired constraints are specified in the "Constraints" block at
the root level of the AMS input file:

   Constraints
      Atom integer
      AtomList integer_list
      FixedRegion string
      Coordinate integer [x|y|z] float?
      Distance (integer){2} float
      All ... [Bonds|Triangles] ...
      Angle (integer){3} float
      Dihedral (integer){4} float
      SumDist (integer){4} float
      DifDist (integer){4} float
      BlockAtoms integer_list
      Block string
      FreezeStrain [xx] [xy] [xz] [yy] [yz] [zz]
      EqualStrain  [xx] [xy] [xz] [yy] [yz] [zz]
   End

"Constraints"
   Type:
      Block

   Description:
      The Constraints block allows geometry optimizations and
      potential energy surface scans with constraints. The constraints
      do not have to be satisfied at the start of the calculation.

The different types constraints are described with examples below.

Note that in principle an arbitrary number of constraints can be
specified and thus combined. However, it is the user's responsibility
to ensure that the specified constraints are actually *compatible with
each other*, meaning that it is theoretically possible to satisfy all
of them at the same time. The AMS driver does **not** detect this kind
of problems, but the optimization will show very unexpected results.
Furthermore, for calculations involving block constraints the
following restrictions apply:

* There should be no other constrained coordinates used together with
  block constraints although this may work in many situation.

* The user should absolutely avoid specifying other constraints that
  include atoms of a frozen block.

Note that a constraint optimization may break symmetry, even if the
constraint(s) themselves would not break symmetry.


Fixed atoms
-----------

The following constraints can be used to fix atoms in space.

"Atom atomIdx"
   Fix the atom with index "atomIdx" at the initial position, as given
   in the "System%Atoms" block.

"AtomList [atomIdx1 .. atomIdxN]"
   Fix all atoms in the list at the initial position, as given in the
   "System%Atoms" block.

"FixedRegion regionName"
   Fix all atoms in a region to their initial positions.

"Coordinate atomIdx [x|y|z] coordValue?"
   Constrain the atom with index "atomIdx" (following the order in the
   "System%Atoms" block) to have a cartesian coordinate ("x", "y" or
   "z") of "coordValue" (given in Angstrom). If the "coordValue" is
   missing, the coordinate will be fixed to its initial value.


Internal degrees of freedom
---------------------------

The following options allow constraints of internal degrees of
freedom, such as distances, angles and dihedral angles.

"Distance atomIdx1 atomIdx2 distValue"
   Constrain the distance between the atoms with index "atomIdx1" and
   "atomIdx2" (following the order in the "System%Atoms" block) to
   "distValue", given in Angstrom.

"Angle atomIdx1 atomIdx2 atomIdx3 angleValue"
   Constrain the angle (1)--(2)--(3) between the atoms with indices
   "atomIdx1-3" (as given by their order in the "System%Atoms" block)
   to "angleValue", given in degrees.

"Dihedral atomIdx1 atomIdx2 atomIdx3 atomIdx4 dihedValue"
   Constrain the dihedral angle (1)--(2)--(3)--(4) between the atoms
   with indices "atomIdx1-4" (as given by their order in the
   "System%Atoms" block) to "dihedValue", given in degrees.

"SumDist atomIdx1 atomIdx2 atomIdx3 atomIdx4 sumDistValue"
   Constrain the sum of the distances R(1,2)+R(3,4) between the atoms
   with indices "atomIdx1-4" (as given by their order in the
   "System%Atoms" block) to "sumDistValue", given in Angstrom.

"DifDist atomIdx1 atomIdx2 atomIdx3 atomIdx4 difDistValue"
   Constrain the difference between the distances R(1,2)-R(3,4) of the
   atoms with indices "atomIdx1-4" (as given by their order in the
   "System%Atoms" block) to "difDistValue", given in Angstrom.

Note that the above constraints do **not** need to be satisfied at the
beginning of the optimization.

The "Distance" keyword described above can be used to set up
constraints for specific pairs of atoms in the system. For some
applications, this input can be cumbersome and one would rather
specify a constraint for **all** bonds of a particular type. Since the
AMS2022.1 release, there is a simplified input for this using the
"All" keyword. The "All" keyword allows generating distance
constraints for matching bonds defined in the System block. The bonds
for which distance constraints will be generated can be selected based
on the elements of the bonded atoms, as well as the bond order defined
in the System block. The following example would constrain all carbon-
carbon single bonds to a distance of 1.4 Angstrom, and all bonds with
hydrogen atoms to their initial distance:

   Constraints
      All single bonds C C to 1.4
      All bonds H *
   End

It is also possible to constrain neighboring bonds and the angle
between them using the combination "All triangles". This is mostly
useful for water molecules, which can then easily be made entirely
rigid:

   All triangles H O H

Warning:

  It is very easy to generate a large number of interdependent
  distance constraints using the "All triangles" keyword.
  Interdependent distance constraints can impair optimizer performance
  and should be avoided. Consider for example the input "All triangles
  H C H" for methane: there are 6 triangles matching this description
  in CH4, generating 18 interdependent distance constraints. A rigid
  CH4 would be better set up using a block constraint.

See below for a full description of the "Constraints%All" keyword.

"Constraints"
   "All"
      Type:
         String

      Recurring:
         True

      Description:
         Fix multiple distances using one the following formats:
         All [bondOrder] bonds at1 at2 [to distance]
         All triangles at1 at2 at3                      The first
         option constrains all bonds between atoms at1 at2 to a
         certain length, while the second - bonds at1-at2 and at2-at3
         as well as the angle between them.                      The
         [bondOrder] can be a number or a string such as single,
         double, triple or aromatic. If it's omitted then any bond
         between specified atoms will be constrained. Atom names are
         case-sensitive and they must be as they are in the Atoms
         block, or an asterisk '*' denoting any atom. If the distance
         is omitted then the bond length from the initial geometry is
         used.                      Important: only the bonds present
         in the system at the start of the simulation can be
         constrained, which means that the bonds may need to be
         specified in the System block.                      Valid
         examples:                      All single bonds C C to 1.4
         All bonds O H to 0.98                      All bonds O H
         All bonds H *                      All triangles H * H


Block constraints
-----------------

Block constraints can be used to treat parts of the System as a rigid
unit in the optimization.

"BlockAtoms [atomIdx1 ... atomIdxN]"
   Creates a block constraint (freezes all internal degrees of
   freedom) for a set of atoms identified by the list of integers
   "[atomIdx1 ...  atomIdxN]". These atom indices refer to the order
   of the atoms in the "System%Atoms" block.

"Block regionName"
   Creates a block constraint (freezes all internal degrees of
   freedom) for a all atoms in a region defined in the "System%Atoms"
   block. Example:

      System
         Atoms
            C  0.0  0.0  0.0    region=myblock
            C  0.0  0.0  1.0    region=myblock
            C  0.0  1.0  0.0
         End
      End
      Constraints
         Block myblock
      End


Lattice constraints
-------------------

For lattice optimizations, the following constraints can be used on
the lattice degrees of freedom:

"FreezeStrain [xx] [xy] [xz] [yy] [yz] [zz]"
   Exclusively for lattice optimizations: Freezes any lattice
   deformation corresponding to a particular component of the strain
   tensor. Accepts a set of strain components [xx, xy, xz, yy, yz, zz]
   to be frozen.

"EqualStrain  [xx] [xy] [xz] [yy] [yz] [zz]"
   Exclusively for lattice optimizations: Accepts a set of strain
   components [xx, xy, xz, yy, yz, zz] which are to be kept equal. The
   applied strain will be determined by the average of the
   corresponding stress tensors components.


Restraints
----------

Not all optimizers support constraints. An alternative is to use so-
called restraints. These are not exact constraints, but rather a
number of springs that pull the system towards the preferred
constraints, see Restraints.


Optimization under pressure / external stress
=============================================

Pressure or non-isotropic external stress can be included in your
simulation via the corresponding engine addons.


Optimization methods
====================

The AMS driver implements a few different geometry optimization
algorithms. It also allows to choose the coordinate space in which the
optimization is performed:

   GeometryOptimization
      Method [Auto | Quasi-Newton | FIRE | L-BFGS | ConjugateGradients | Dimer]
      CoordinateType [Auto | Delocalized | Cartesian]
   End

"GeometryOptimization"
   "Method"
      Type:
         Multiple Choice

      Default value:
         Auto

      Options:
         [Auto, Quasi-Newton, FIRE, L-BFGS, ConjugateGradients, Dimer]

      GUI name:
         Optimization method

      Description:
         Select the optimization algorithm employed for the geometry
         relaxation. Currently supported are:
         the Hessian-based Quasi-Newton-type BFGS algorithm,
         the fast inertial relaxation method (FIRE),
         the limited-memory BFGS method,                       and the
         conjugate gradients method. The default is to choose an
         appropriate method automatically based on the engine's speed,
         the system size and the supported optimization options.

   "CoordinateType"
      Type:
         Multiple Choice

      Default value:
         Auto

      Options:
         [Auto, Delocalized, Cartesian]

      GUI name:
         Optimization space

      Description:
         Select the type of coordinates in which to perform the
         optimization. 'Auto' automatically selects the most
         appropriate CoordinateType for a given Method.
         If 'Auto' is selected, Delocalized coordinates will be used
         for the Quasi-Newton method, while Cartesian coordinates will
         be used for all other methods.

We **strongly** advise leaving both the "Method" as well as the
"Coordinate" type on the "Auto" setting. There are many restrictions
as to which optimizer and coordinate type can be used together with
which kind of optimization. The following (roughly) sketches the
compatibility of the different optimizers and options:

+-----------------------+----------------------------------------+----------------+------------------------+
| Optimizer             | Constraints                            | Lattice opt.   | Coordinate types       |
|=======================|========================================|================|========================|
| Quasi-Newton          | All, except strain                     | Yes            | All                    |
+-----------------------+----------------------------------------+----------------+------------------------+
| FIRE                  | Fixed coordinates, distances, strain   | Yes            | Cartesian              |
+-----------------------+----------------------------------------+----------------+------------------------+
| L-BFGS                | No                                     | Yes            | Cartesian              |
+-----------------------+----------------------------------------+----------------+------------------------+
| Conjugate gradients   | No                                     | No             | Cartesian              |
+-----------------------+----------------------------------------+----------------+------------------------+

Furthermore for optimal performance the optimizer should be chosen
with the speed of the engine: a faster engine in combination should
use an optimizer with little overhead (FIRE), while slower engines
should use optimizers that strictly minimize the number of steps
(Quasi-Newton). This is all handled automatically by default, and we
recommend changing "Method" and "Coordinate" only in case there are
problems with the automatic choice.

The following subsections list the strengths and weaknesses of the
individual optimizers in some more detail, motivating why which
optimizer is chosen automatically under which circumstances.


Quasi-Newton
------------

This optimizer implements a quasi Newton approach [1] [2] [3], using
the Hessian for computing changes in the geometry so as to reach a
local minimum. The Hessian itself is typically approximated in the
beginning and updated in the process of optimization. It uses
delocalized coordinates by default both for molecules and periodic
systems. The molecular part is based mainly on the work by Marcel
Swart [4]. Cartesian coordinates are used in the presence of an
external electric field and/or frozen atom constraints.

The Quasi-Newton (QN) optimizer supports all types of constraints and
can be used for both molecular and periodic systems, including lattice
optimizations. For cases where the optimization can be performed in
delocalized coordinates, the number of steps taken to reach the local
minimum is usually smaller than when optimizing in Cartesian ones.
For fast compute engines, the overhead of the QN optimizer can become
a bottleneck of the calculation, thus a more light-weight optimizer
such as FIRE may give an better overall performance. In principle, a
QN optimization in delocalized coordinates may run out of memory for a
very large system (say over 1000 atoms) because of the SVD step.
However, since it is going to be used for a moderate-to-slow engine we
still recommend sticking to it for the benefit of fewer steps. Because
of these properties the QN optimizer is the default in AMS for all
kinds of optimizations with moderate and slow engines, such DFTB and
ADF.  It is also used as the optimizer back-end for the PES scan task,
the transition state search as well as the calculation of the elastic
tensor.

Details of the Quasi-Newton optimizer are configured in a dedicated
block:

   GeometryOptimization
      Quasi-Newton
         MaxGDIISVectors integer
         Step
            TrustRadius float
            VaryTrustRadius Yes/No
         End
         UpdateTSVectorEveryStep Yes/No
      End
   End

"GeometryOptimization"
   "Quasi-Newton"
      Type:
         Block

      Description:
         Configures details of the Quasi-Newton geometry optimizer.

      "MaxGDIISVectors"
         Type:
            Integer

         Default value:
            0

         Description:
            Sets the maximum number of GDIIS vectors. Setting this to
            a number >0 enables the GDIIS method.

      "Step"
         Type:
            Block

         Description:
         "TrustRadius"
            Type:
               Float

            Description:
               Initial value of the trust radius.

         "VaryTrustRadius"
            Type:
               Bool

            Description:
               Whether to allow the trust radius to change during
               optimization. By default True during energy
               minimization and False during transition state search.

      "UpdateTSVectorEveryStep"
         Type:
            Bool

         Default value:
            Yes

         GUI name:
            Update TSRC vector every step

         Description:
            Whether to update the TS reaction coordinate at each step
            with the current eigenvector.

The Quasi-Newton optimizer uses the Hessian to compute the step of the
geometry optimization. The Hessian is typically approximated in the
beginning and then updated during the optimization. A very good
initial Hessian can therefore increase the performance of the
optimizer and lead to faster and more stable convergence. The choice
of the initial Hessian can be configured in a dedicated block:

   GeometryOptimization
      InitialHessian
         File string
         Type [Auto | UnitMatrix | Swart | FromFile | Calculate | CalculateWithFastEngine]
      End
   End

"GeometryOptimization"
   "InitialHessian"
      Type:
         Block

      Description:
         Options for initial model Hessian when optimizing systems
         with the Quasi-Newton method.

      "File"
         Type:
            String

         GUI name:
            Initial Hessian from

         Description:
            KF file containing the initial Hessian (or the results
            dir. containing it). This can be used to load a Hessian
            calculated in a previously with the [Properties%Hessian]
            keyword.

      "Type"
         Type:
            Multiple Choice

         Default value:
            Auto

         Options:
            [Auto, UnitMatrix, Swart, FromFile, Calculate,
            CalculateWithFastEngine]

         GUI name:
            Initial Hessian

         Description:
            Select the type of initial Hessian. Auto: let the program
            pick an initial model Hessian. UnitMatrix: simplest
            initial model Hessian, just a unit matrix in the
            optimization coordinates. Swart: model Hessian from M.
            Swart. FromFile: load the Hessian from the results of a
            previous calculation (see InitialHessian%File). Calculate:
            compute the initial Hessian (this may be computationally
            expensive and it is mostly recommended for
            TransitionStateSearch calculations).
            CalculateWithFastEngine: compute the initial Hessian with
            a faster engine.

While there are some options for the construction of approximate model
Hessians, the best initial Hessians are often those calculated
explicitly at a lower level of theory, e.g. the real DFTB Hessian can
be used the initial Hessian for an optimization with the more accurate
BAND engine. Using the "CalculateWithFasterEngine" keyword can be used
to automatically chose a fast engine at a lower level of theory.  What
the lower level of theory is depends on the main engine used in the
calculation: DFTB with the GFN1-xTB model is used as the lower level
of theory for relatively slow engines, e.g. DFT based engines. For
semi-empirical engines like DFTB or MOPAC, the lower level of theory
is currently UFF. If more control over the lower level engine is
needed, the initial Hessian can be calculated with a user defined
engine and then loaded from file, see this example.


FIRE
----

The Fast Inertial Relaxation Engine [5] based optimizer has basically
no overhead per step, so that the speed of the optimization purely
depends on the performance of the used compute engine. As such it is a
good option for large systems or fast compute engines, where the
overhead of the Quasi-Newton optimizer would be significant. Note that
is also supports fixed atom constraints and coordinate constraints (as
long as the value of the constrained coordinate is already satisfied
in the input geometry), distance constraints, as well as lattice
optimizations (with strain constraints).

FIRE is selected as the default optimizer for fast compute engines if
it is compatible with all other settings of the optimization (i.e. no
unsupported constraints or coordinate types).

Note:

  FIRE is a very robust optimizer. In case of convergence problems
  with the other methods, it is a good idea to see if the optimization
  converges with FIRE. If it does not, it is very likely that the
  problem is not the optimizer but the shape of the potential energy
  surface ...

The details of the FIRE optimizer are configured in a dedicated block.
It is quite easy to make the optimization numerically unstable when
tweaking these settings, so we strongly recommend leaving everything
at the default values.

   GeometryOptimization
      FIRE
         AllowOverallRotation Yes/No
         AllowOverallTranslation Yes/No
         MapAtomsToUnitCell Yes/No
         NMin integer
         alphaStart float
         dtMax float
         dtStart float
         fAlpha float
         fDec float
         fInc float
         strainMass float
      End
   End

"GeometryOptimization"
   "FIRE"
      Type:
         Block

      Description:
         This block configures the details of the FIRE optimizer. The
         keywords name correspond the the symbols used in the article
         describing the method, see PRL 97, 170201 (2006).

      "AllowOverallRotation"
         Type:
            Bool

         Default value:
            Yes

         Description:
            Whether or not the system is allowed to freely rotate
            during the optimization. This is relevant when optimizing
            structures in the presence of external fields.

      "AllowOverallTranslation"
         Type:
            Bool

         Default value:
            No

         Description:
            Whether or not the system is allowed to translate during
            the optimization. This is relevant when optimizing
            structures in the presence of external fields.

      "MapAtomsToUnitCell"
         Type:
            Bool

         Default value:
            No

         Description:
            Map the atoms to the central cell at each geometry step.

      "NMin"
         Type:
            Integer

         Default value:
            5

         Description:
            Number of steps after stopping before increasing the time
            step again.

      "alphaStart"
         Type:
            Float

         Default value:
            0.1

         Description:
            Steering coefficient.

      "dtMax"
         Type:
            Float

         Default value:
            1.0

         Unit:
            Femtoseconds

         Description:
            Maximum time step used for the integration. For ReaxFF and
            APPLE&P, this value is reduced by 50%.

      "dtStart"
         Type:
            Float

         Default value:
            0.25

         Unit:
            Femtoseconds

         Description:
            Initial time step for the integration.

      "fAlpha"
         Type:
            Float

         Default value:
            0.99

         Description:
            Reduction factor for the steering coefficient.

      "fDec"
         Type:
            Float

         Default value:
            0.5

         Description:
            Reduction factor for reducing the time step in case of
            uphill movement.

      "fInc"
         Type:
            Float

         Default value:
            1.1

         Description:
            Growth factor for the integration time step.

      "strainMass"
         Type:
            Float

         Default value:
            0.5

         Description:
            Fictitious relative mass of the lattice degrees of
            freedom. This controls the stiffness of the lattice
            degrees of freedom relative to the atomic degrees of
            freedom, with smaller values resulting in a more
            aggressive optimization of the lattice.


Limited-memory BFGS
-------------------

AMS also offers an L-BFGS based geometry optimizer. It usually
converges faster than FIRE, but does not support constrained
optimizations. For periodic systems it can be quite good for lattice
optimizations. The new implementation has not been thoroughly tested
yet, therefore never selected automatically. For large systems and
fast engines you may want to disable symmetry: simply the detection of
(non-existing) symmetry may be a huge overhead.

   GeometryOptimization
      HessianFree
         Step
            MaxCartesianStep float
            MinRadius float
            TrialStep float
            TrustRadius float
         End
      End
   End

"GeometryOptimization"
   "HessianFree"
      Type:
         Block

      Description:
         Configures details of the Hessian-free (conjugate gradients
         or L-BFGS) geometry optimizer.

      "Step"
         Type:
            Block

         Description:
         "MaxCartesianStep"
            Type:
               Float

            Default value:
               0.1

            Unit:
               Angstrom

            Description:
               Limit on a single Cartesian component of the step.

         "MinRadius"
            Type:
               Float

            Default value:
               0.0

            Unit:
               Angstrom

            Description:
               Minimum value for the trust radius.

         "TrialStep"
            Type:
               Float

            Default value:
               0.0005

            Unit:
               Angstrom

            Description:
               Length of the finite-difference step when determining
               curvature. Should be smaller than the step convergence
               criterion.

         "TrustRadius"
            Type:
               Float

            Default value:
               0.2

            Unit:
               Angstrom

            Description:
               Initial value of the trust radius.


Conjugate gradients
-------------------

AMS also offers a conjugate gradients based geometry optimizer, as it
was also implemented in the pre-2018 releases of the DFTB program.
However, it is usually slightly slower than FIRE, and supports neither
lattice nor constrained optimizations. It is therefore never selected
automatically, and we do not recommend using it. Like L-BFGS, the
conjugate gradients optimizer is also configured in the "HessianFree"
block, see L-BFGS section above for details.


Dimer method for TS seach
-------------------------

The dimer method [6] works differently from the quasi-Newton. It is
designed to always follow the lowest mode from the initial geometry
without requiring the Hessian. The mode (a.k.a. the dimer vector) is
defined by two close-lying points (midpoint and endpoint) on the
potential energy surface. The vector's direction is found by
minimizing the PES curvature (finite-difference second derivative of
the energy along the vector) on the hypersphere around the midpoint.
This is achieved by iteratively rotating the vector until the
estimated rotation angle drops below the "AngleThreshold" value. These
are called rotation iterations. After that, a translation iteration is
performed that maximized the energy along the dimer vector and
minimizes it in all other directions. The optimization is considered
converged when the largest energy gradient component drops below the
"Convergence%Gradient" value.

Both the rotation and the translation are performed in Cartesian
coordinates using the L-BFGS method as described in [7]. The L-BFGS
trust radii are set by the "RotationTrustRadius" and
"TranslationTrustRadius" keywords, respectively. The TS search can be
restricted to a subset of atomic coordinates defined by the "Region"
keyword, in which case only atoms of the selected region are moved
during rotations.

Each translation iteration requires two engine evaluations: one for
the midpoint and one for the endpoint. Each rotation iteration
additionally requires one or two engine evaluations depending on the
"ExtrapolateForces" setting. It should be noted that with the default
settings the rotations are usually necessary only at the beginning of
the optimization. At the later steps the dimer vector (the direction
of the search) does not usually change much.

By default, the initial value of the search vector is chosen randomly.
This can be changed using either a "ReactionCoordinate" or a "System
final" input block described in the Transition state search section.

   GeometryOptimization
      Dimer
         AngleThreshold float
         DimerDelta float
         ExtrapolateForces Yes/No
         LBFGSMaxVectors integer
         MaxRotationIterations integer
         Region string
         RotationTrustRadius float
         TranslationTrustRadius float
      End
   End

"GeometryOptimization"
   "Dimer"
      Type:
         Block

      Description:
         Options for the Dimer method for transition state search.

      "AngleThreshold"
         Type:
            Float

         Default value:
            1.0

         Unit:
            Degree

         Description:
            The rotation is considered converged when the the rotation
            angle falls below the specified threshold.

      "DimerDelta"
         Type:
            Float

         Default value:
            0.01

         Unit:
            Angstrom

         Description:
            Euclidian distance between the midpoint and the endpoint.

      "ExtrapolateForces"
         Type:
            Bool

         Default value:
            Yes

         Description:
            Set to false to call engine to calculate forces at the
            extrapolated rotation angle instead of extrapolating them.

      "LBFGSMaxVectors"
         Type:
            Integer

         Default value:
            10

         Description:
            Max number of vectors for the L-BFGS algorithm to save.

      "MaxRotationIterations"
         Type:
            Integer

         Default value:
            10

         Description:
            Maximum number of rotation iterations for a single
            translation step.

      "Region"
         Type:
            String

         Default value:
            *

         Description:
            Include only atoms of the specified region(s) in the
            rotations, which allows searching for a transition state
            involving selected atoms only.

      "RotationTrustRadius"
         Type:
            Float

         Default value:
            0.1

         Unit:
            Angstrom

         Description:
            L-BFGS trust radius during rotation iterations.

      "TranslationTrustRadius"
         Type:
            Float

         Default value:
            0.1

         Unit:
            Angstrom

         Description:
            L-BFGS trust radius during translation iterations.


Troubleshooting
===============


Failure to converge
-------------------

First of all one should look how the energy changed during the latest
ten or so iterations. If the energy is decreasing more or less
consistently, possibly with occasional jumps, then there is probably
nothing wrong with the optimization. This behavior is typical in the
cases when the starting geometry was far away from the minimum and the
optimization has a long way to go. Just increase the allowed number of
iterations, restart from the latest geometry and see if the
optimization converges.

If the energy oscillates around some value and the energy gradient
hardly changes then you may need to look at the calculation setup.

The success of geometry optimization depends on the accuracy of the
calculated forces. The default accuracy settings are sufficient in
most cases. There are, however, cases when one has to increase the
accuracy in order to get geometry optimization converged. First of
all, this may be necessary if you tighten the optimization convergence
criteria. In some cases it may be necessary to increase the accuracy
also for the default criteria. Please refer to the engine manuals for
instructions on how to increase the accuracy of an engine's energies
and gradients. Often this is done with the "NumericalQuality" keyword
in the engine input.

A geometry optimization can also fail to converge because the
underlying potential energy surface is problematic, e.g. it might be
discontinuous or not have a minimum at which the gradients vanish.
This often indicates real problems in the calculation setup, e.g. an
electronic structure that changes fundamentally between subsequent
steps in the optimization. In these cases it is advisable to run a
single point calculation at the problematic geometries and carefully
check if the results are physically actually sensible.

Finally it can also be a technical problem with the specific
optimization method used. In these cases switching to another method
could help with convergence problems. We recommend first trying the
FIRE optimizer, as it is internally relatively simple and stable.


Converged, but did not find a minimum on the PES
------------------------------------------------

If Normal Modes or  PES point characterization were requested, the PES
character of the final optimized geometry is checked.

While you would generally expect the final optimized geometry to be a
local minimum (i.e. no imaginary frequencies), it can happen for the
optimization procedure to converge to a saddle point instead (i.e. a
stationary point on the PES with one or more imaginary frequency). In
that case you will get the following error message:

   ERROR: Geometry optimization failed! (Converged, but did not find a minimum on the PES.)

Depending on your application, (small) imaginary frequencies might not
be a problem. But if you want to get rid of all imaginary frequencies,
these are a few recommendations:

* Try using a tighter gradient convergence threshold
  ("GeometryOptimization%Convergence%Gradient"). E.g. 1.0E-4 (instead
  of the default 1.0E-3).

* Try using the MaxRestart option. Reasonable values for MaxRestart
  range from 2 to 20 (if there are many imaginary frequencies, a
  larger value might be needed). When using the "MaxRestart" option,
  you might try increasing "RestartDisplacement".

* Open the results with the **AMSSpectra** GUI module and visually
  inspect the negative/imaginary vibrational modes. You can then
  manually perturb the input geometry along those modes and re-run the
  optimization.

* If the system is symmetric and symmetry is enabled, the optimizer
  will try not to break the symmetry of the system. It might be that
  in order to reach the local minimum, symmetry must be broken. Try
  disabling symmetry ("UseSymmetry False") and slightly perturb the
  atomic coordinates in order to break the symmetry (or use the
  PerturbCoordinates option.)


Restarting a geometry optimization
----------------------------------

During a running optimization the system's geometry is written out to
the AMS driver's output file "ams.rkf" after every step (in the
"Molecule" section). This means that crashed or otherwise canceled
geometry optimizations can be restarted by simply loading the last
frame from there using the "LoadSystem" keyword, see its documentation
in the system definition section of this manual:

   LoadSystem File=my_crashed_GO.results/ams.rkf

This can of course also be used to continue an optimization but e.g.
with tighter convergence criteria or a different optimizer, as it
essentially starts a new geometry optimization from the previous
geometry, and does not propagate any information internal to the
optimizer (e.g.  the approximate Hessian for the Quasi-Newton
optimizer or the FIRE velocities) to the new job. As such it might
take a few more steps to convergence than if the original job had
continued, but allows for additional flexibility.


Restart an Engine from file
~~~~~~~~~~~~~~~~~~~~~~~~~~~

   EngineRestart string

"EngineRestart"
   Type:
      String

   Description:
      The path to the file from which to restart the engine.
      Should be a proper engine result file (like adf.rkf, band.rkf
      etc), or the name of the results directory containing it.

EngineRestart should typically be combined with LoadSystem and
depending on whether one wants to change the engine configuration may
be combined LoadEngine. For example, in case of ADF a restart of a
geometry optimization (use TAPE13 instead of adf.rkf in case of a
crash) may look like:

   Task GeometryOptimization
   LoadSystem
       File crashed_or_not_crashed.results/ams.rkf
   End
   EngineRestart not_crashed.results/adf.rkf
   LoadEngine not_crashed.results/adf.rkf

[1] L. Versluis and T. Ziegler, *The determination of Molecular
    Structure by Density Functional Theory*, Journal of Chemical
    Physics 88, 322 (1988)

[2] L. Versluis, The determination of molecular structures by the HFS
    method, PhD thesis, University of Calgary, 1989

[3] L. Fan and T. Ziegler, *Optimization of molecular structures by
    self consistent and non-local density functional theory*, Journal
    of Chemical Physics 95, 7401 (1991)

[4] M. Swart and F.M. Bickelhaupt, *Optimization of strong and weak
    coordinates*, International Journal of Quantum Chemistry 106, 2536
    (2006)

[5] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch,
    *Structural Relaxation Made Simple*, Physical Review Letters 97,
    170201 (2006)

[6] 7. Henkelman and H. Jonsson, *A dimer method for finding saddle
       points on high dimensional potential surfaces using only first
       derivatives*, Journal of Chemical Physics 111, 7010 (1999)

[7] 10. Kästner and P. Sherwood, *Superlinearly converging dimer
        method for transition state search*, Journal of Chemical
        Physics 128, 014106 (2008)
