Geometry optimization¶

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

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.

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

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.

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.

EngineRestart string

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