# 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
Step float
End
MaxIterations integer
CalcPropertiesOnlyIfConverged [True | False]
OptimizeLattice [True | False]
Pressure float
KeepIntermediateResults [True | False]
End

GeometryOptimization
Type: Block Configures details of the geometry optimization and transition state searches.
Convergence
Type: Block Convergence is monitored for two items: the energy and the Cartesian gradients. Convergence criteria can be specified separately for each of these items.
Energy
Type: Float 1e-05 Hartree The criterion for changes in the energy.
Gradients
Type: Float 0.001 Hartree/Angstrom The criterion for changes in the gradients.
Step
Type: Float 0.001 Angstrom The maximum Cartesian step allowed for a converged geometry.

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

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 The maximum number of geometry iterations allowed to converge to the desired structure.
CalcPropertiesOnlyIfConverged
Type: Bool True 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.

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.

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

GeometryOptimization
OptimizeLattice
Type: Bool False Whether to also optimize the lattice for periodic structures. This is currently only supported with the Quasi-Newton and SCMGO optimizers.
Pressure
Type: Float 0.0 Optimize the structure under pressure (this will only have an effect if you are optimizing the lattice vectors). Currently only working in combination with the Quasi-Newton optimizer. For phase transitions you may consider disabling or breaking the symmetry.

Finally the GeometryOptimization block also contains some technical options:

GeometryOptimization
KeepIntermediateResults
Type: Bool False 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
Coordinate integer [x|y|z] float?
Distance (integer){2} float
Angle (integer){3} float
Dihedral (integer){4} float
BlockAtoms integer_list
Block string
End

Atom atomIdx
Fix the atom with index atomIdx at the initial position, as given in the System%Atoms block.
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.
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.
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 blockName

Creates a block constraint (freezes all internal degrees of freedom) for a set of atoms identified by a tagging string blockName in the System%Atoms block. The tag is attached to element symbol and separated by a dot. Example:

System
Atoms
C.myblock  0.0  0.0  0.0
C.myblock  0.0  0.0  1.0
C          0.0  1.0  0.0
End
End
Constraints
Block myblock
End


Note that the coordinate, distance, angle, and dihedral constraints do not need to be satisfied at the beginning of the optimization.

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 these kinds 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.

## 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 | SCMGO | FIRE | ConjugateGradients]
CoordinateType [Auto | Delocalized | Cartesian]
End

GeometryOptimization
Method
Type: Multiple Choice Auto [Auto, Quasi-Newton, SCMGO, FIRE, ConjugateGradients] Select the optimization algorithm employed for the geometry relaxation. Currently supported are: the Hessian-based Quasi-Newton-type BFGS algorithm, the experimental SCMGO optimizer, the fast inertial relaxation method (FIRE), 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 Auto [Auto, Delocalized, Cartesian] Select the type of coordinates in which to perform the optimization. If ‘Auto’, delocalized coordinates will be used for molecular systems, while cartesian coordinates will be used for periodic systems. Optimization in delocalized coordinates [Delocalized] can only be used for geometry optimizations or transition state searches of molecular systems with the Quasi-Newton method. The experimental SCMGO optimizer supports [Delocalized] coordinates for both molecular and periodic systems.

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 compatibilty of the different optimizers and options:

Optimizer Constraints Lattice opt. Coordinate types
Quasi-Newton Yes Yes All (molecules)
Cartesian (periodic)
SCMGO No Yes Delocalized
FIRE Fixed atoms and coordinates Yes Cartesian
Conjugate gradients No No Cartesian

Furthermore for optimal performance the optimizer should be chosen with the speed of the engine in mind: Faster engine should use an optimizer with little overhead (FIRE), while slower engines should use optimizers that strictly minimize the number of steps (Quasi-Newton, SCMGO). 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-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. For molecules it uses delocalized coordinates by default, based mainly on the work by Marcel Swart . For periodic systems the optimization is performed in cartesian coordinates instead.

The Quasi-Newton optimizer supports all types of constraints and can be used for both molecular and periodic systems, including lattice optimizations. For molecular systems, where the optimization can be performed in delocalized coordinates, the number of steps taken to reach the local minimum is quite small. For large systems (on the order of hundreds of atoms) and fast compute engines, the overhead of the Quasi-Newton optimizer is likely the bottleneck of the calculation, and more light-weight optimizers like FIRE will give an overall better performance. We do not recommend using the Quasi-Newton optimizer for systems >1000 atoms. Because of these properties the Quasi-Newton optimizer is the default in AMS for all kinds of optimizations of small and medium sized systems. It is also used as the optimizer backend 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
End
End
End

GeometryOptimization
Quasi-Newton
Type: Block Configures details of the Quasi-Newton geometry optimizer.
MaxGDIISVectors
Type: Integer 0 Sets the maximum number of GDIIS vectors. Setting this to a number >0 enables the GDIIS method.
Step
Type: Block
TrustRadius
Type: Float 0.2 Initial value of the trust radius.

The Quasi-Newton optimizer uses the Hessian to compute the next 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]
End
End

GeometryOptimization
InitialHessian
Type: Block Options for initial model Hessian when optimizing systems with either the Quasi-Newton or the SCMGO method.
File
Type: String KF file containing the initial Hessian. This can be used to load a Hessian calculated in a previously with the [Properties%Hessian] keyword.
Type
Type: Multiple Choice Auto [Auto, UnitMatrix, Swart, FromFile] Selects the type of the initial model Hessian, or load the Hessian from the results of a previous calulation.

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, see this example.

### SCMGO¶

The SCMGO optimizer is a new implementation of a Quasi-Newton style optimizer working in delocalized coordinates. In the 2018 release of the Amsterdam Modeling Suite it is still considered experimental and therefore never selected automatically. However, for molecules and fully connected periodic systems it already shows a quite good performance, and could be a reasonable alternative to the classic Quasi-Newton optimizer, which can not use the more efficient delocalized coordinates for periodic systems.

GeometryOptimization
SCMGO
ContractPrimitives [True | False]
NumericalBMatrix [True | False]
Step
VariableTrustRadius [True | False]
End
logSCMGO [True | False]
testSCMGO [True | False]
End
End

GeometryOptimization
SCMGO
Type: Block Configures details SCMGO.
ContractPrimitives
Type: Bool True Form non-redundant linear combinations of primitive coordinates sharing the same central atom
NumericalBMatrix
Type: Bool False Calculation of the B-matrix, i.e. Jacobian of internal coordinates in terms of numerical differentiations
Step
Type: Block
TrustRadius
Type: Float 0.2 Initial value of the trust radius.
VariableTrustRadius
Type: Bool True Whether or not the trust radius can be updated during the optimization.
logSCMGO
Type: Bool False Verbose output of SCMGO internal data
testSCMGO
Type: Bool False Run SCMGO in test mode.

Note that SCMGO also supports different initial Hessians, and uses the same options for the initial Hessian as the Quasi-Newton optimizer, see above.

### FIRE¶

The Fast Inertial Relaxation Engine  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), as well as lattice optimizations.

FIRE is the default optimizer for systems >1000 atoms. For smaller systems and fast compute engines it is selected 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
NMin integer
alphaStart float
dtMax float
dtStart float
fAlpha float
fDec float
fInc float
strainMass float
End
End

GeometryOptimization
FIRE
Type: Block 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).
NMin
Type: Integer 5 Number of steps after stopping before increasing the time step again.
alphaStart
Type: Float 0.1 Steering coefficient.
dtMax
Type: Float 1.25 Femtoseconds Maximum time step used for the integration.
dtStart
Type: Float 0.25 Femtoseconds Initial time step for the integration.
fAlpha
Type: Float 0.99 Reduction factor for the steering coefficient.
fDec
Type: Float 0.5 Reduction factor for reducing the time step in case of uphill movement.
fInc
Type: Float 1.1 Growth factor for the integration time step.
strainMass
Type: Float 0.5 Ficticious 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.

Note that neither the energy change per step, nor the step size are reliable convergence criteria for the FIRE optimizer. Only the gradient convergence criterium (set with the Converge%Gradients key) is used by FIRE to determine when the optimization has converged.

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 less robust, and has a larger overhead than FIRE, and supports neither lattice nor constrained optimizations. It is therefore never selected automatically, and we do not recommend using it.

GeometryOptimization
Step
End
End
End

GeometryOptimization
ConjugateGradients
Type: Block Configures details of the conjugate gradients geometry optimizer.
Step
Type: Block
MinRadius
Type: Float 0.0
TrustRadius
Type: Float 0.2 Initial value of the trust radius.

## 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, stable and has the nice property that the energy can never increase during an optimization.

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