Geometry optimization¶

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 [4]. 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
         TrustRadius float
      End
   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
Default value:	0.2
Description:	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
Description:	Options for initial model Hessian when optimizing systems with either the Quasi-Newton or the SCMGO method.

File

Type:	String
Description:	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
Default value:	Auto
Options:	[Auto, UnitMatrix, Swart, FromFile]
Description:	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
         TrustRadius float
         VariableTrustRadius [True | False]
      End
      logSCMGO [True | False]
      testSCMGO [True | False]
   End
End

GeometryOptimization

SCMGO

Type:	Block
Description:	Configures details SCMGO.

ContractPrimitives

Type:	Bool
Default value:	True
Description:	Form non-redundant linear combinations of primitive coordinates sharing the same central atom

NumericalBMatrix

Type:	Bool
Default value:	False
Description:	Calculation of the B-matrix, i.e. Jacobian of internal coordinates in terms of numerical differentiations

Step

Type:	Block
Description:

TrustRadius

Type:	Float
Default value:	0.2
Description:	Initial value of the trust radius.

VariableTrustRadius

Type:	Bool
Default value:	True
Description:	Whether or not the trust radius can be updated during the optimization.

logSCMGO

Type:	Bool
Default value:	False
Description:	Verbose output of SCMGO internal data

testSCMGO

Type:	Bool
Default value:	False
Description:	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 [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), 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).

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

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.25
Unit:	Femtoseconds
Description:	Maximum time step used for the integration.

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

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 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
   ConjugateGradients
      Step
         MinRadius float
         TrustRadius float
      End
   End
End

GeometryOptimization

ConjugateGradients

Type:	Block
Description:	Configures details of the conjugate gradients geometry optimizer.

Step

Type:	Block
Description:

MinRadius

Type:	Float
Default value:	0.0
Description:

TrustRadius

Type:	Float
Default value:	0.2
Description:	Initial value of the trust radius.

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.