# 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 systems. 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.

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:

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.The maximum Cartesian nuclear gradient is smaller than

`Convergence%Gradient`

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

`Convergence%Gradient`

.The maximum Cartesian step is smaller than

`Convergence%Step`

.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 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 to perform constrained optimizations, where a number of specified degrees of freedom are fixed to particular values.

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.

### 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.

References

- 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
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
Kästner and P. Sherwood,

*Superlinearly converging dimer method for transition state search*, Journal of Chemical Physics 128, 014106 (2008)