Geometry Optimization

The Geometry Optimizations procedure in ADF is based on 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 initialized and updated in the process of optimization. The default optimization strategy uses delocalized coordinates and it is mainly based on the work by Marcel Swart [4].

See also

In case of geometry convergence problems, see Geometry Optimization troubleshooting

You can perform a Geometry Optimization using default settings by simply specifying the GEOMETRY key block:

GEOMETRY
End

Several subkeys in the GEOMETRY block can be used for control of the Geometry Optimization procedure and related parameters.

GEOMETRY
   Converge {E=TolE} {Grad=TolG} {Rad=TolR} {Angle=TolA}
   Iterations Niter
   Inithess inithessian.file
   Hessupd HessUpdate
   Optim {Delocal / Cartesian / Internal} {All / Selected}
   GEigenvalueThreshold threshold
   Branch {New / Old}
   Step {Rad=MaxRadStep} {Angle=MaxAngleStep} {TrustRadius=MaxRadius}
   DIIS {N=NVect} {CYC=Ncyc}
   Externprogram externprog.exe coord=coords.inp energy=energy.out grad=grads.out
End

Convergence

A geometry optimization is considered converged when all the following criteria are met (this information is also printed to the logfile):

  1. the difference between the bond energy at the current geometry and at the previous geometry step is smaller than \(\text{TolE}\)
  2. the maximum Cartesian nuclear gradient is smaller than \(\text{TolG}\)
  3. the root mean square (RMS) of the Cartesian nuclear gradients is smaller than \(\frac{2}{3} \text{TolG}\)
  4. the maximum Cartesian step is smaller than \(\text{TolR}\)
  5. the root mean square (RMS) of the Cartesian steps is smaller than \(\frac{2}{3} \text{TolR}\)

note: if the maximum and RMS gradients are 10 times smaller than the convergence criterion, then criteria 4 and 5 are ignored

The convergence criteria can be changed via the Converge key:

GEOMETRY
   Converge {E=TolE} {Grad=TolG} {Rad=TolR} {Angle=TolA}
END
TolE
The criterion for changes in the energy, in Hartree. Default: 1e-3.
TolG
The criterion for the nuclear gradients, in Hartree/Angstrom. Default: 1e-3.
TolR
Refers to changes in the Cartesian coordinates or bond lengths, depending on in what coordinates you optimize, in angstrom. Default: 1e-2.
TolA
Refers to changes in bond- and dihedral angles, in degrees. This is only meaningful if optimization takes place in Z-matrix coordinates. Default: 0.5 degree.

Remarks:

  • 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, a loose threshold may yield geometries that are far from the minimum as regards 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 ADF are intended to be reasonable for most applications, but inevitably situations may arise where they are inadequate.
  • The technical numerical accuracy of the calculation (e.g. numerical integration) should somehow match the required level of convergence in gradients (strict convergence criteria may require high numerical accuracy). A simple way of changing the numerical accuracy is through the NumericalQuality keyword.
  • The convergence threshold for the coordinates (TolL, TolA) 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 the program-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.

Iterations

GEOMETRY
   Iterations Niter
END
Iterations
Niter
The maximum number of geometry iterations allowed to locate the desired structure. The default is max(30,2*Nfree), where Nfree is the number of free variables, which is typically close to 3*N(atoms). This is a fairly large number. If the geometry has not converged (at least to a reasonable extent) within that many iterations you should sit down and consider the underlying cause rather than simply increase the allowed number of cycles and try again.

Optimization strategy

Optim

GEOMETRY
   Optim {Delocal / Cartesian / Internal} {All / Selected}
   GEigenvalueThreshold threshold
END
Optim
Delocal
Optimization in delocalized coordinates [4] (Delocal) is the default in geometry optimizations, transition state searches, and (linear) transits. You can use constraints (see the CONSTRAINTS key) and restraints, and you can supply the atoms in Z-matrix coordinates. You can not use dummy atoms, ghost atoms, or alternative elements.
Cartesian / Internal

Cartesian or Zmatrix (equivalently: internal) specifies the type of coordinates in which the minimization is carried out. By default the coordinate type is applied that was used in the ATOMS key for the input of the (initial) atomic positions. (Cartesian if atoms were input in zcart format). Note that Zmatrix optimizations can only be done with the old branch of optimizations.

Cartesian optimization is allowed if the atoms were input in Z-matrix format, but no constraints (see the key GEOVAR) can then be used: all coordinates are optimized. An attempt to explicitly freeze variables may result in an error abort. Optimization in Z-matrix coordinates is not allowed if only Cartesian were supplied in atoms: the program does not construct a Z-matrix by itself. One should then use the zcart format: give Cartesian coordinates and supply the structure of the Z-matrix. Again, in this case you cannot use constraints.

Selected
For use with old branch of optimizations. Only those coordinates are optimized that are defined with the key GEOVAR.

All (The default value) means that in principle all atomic coordinates will be varied. With the new branch of optimization the key CONSTRAINTS one may set constraints on the optimization. With the old branch of optimization one can use the key GEOVAR to set constraints.

GEigenvalueThreshold
threshold
When performing optimization in delocalized coordinates a singular value decomposition of the Wilson B-matrix is obtained. Left-eigenvectors of the B-matrix corresponding to non-zero singular values are used to determine the delocalized coordinates. ‘Non-zero’ in this case means larger than sqrt(threshold). The square root is taken because left-eigenvectors of the B-matrix correspond to eigenvectors of the G = B Bt symmetric matrix and the singular values of B are related to eigenvalues of G as \(\lambda_i(G) = {\sigma_i(B)}^2\). The default value is 10-15. Note: in rare cases, ADF may generate a “Suspicious delocalized coordinates” warning. In such a case, check the output file for more details. For non-linear molecules there should be no more than 3N-6 eigenvalues above the threshold. There may be fewer than that if block constraints are used. If the number of delocalized coordinates is higher than 3N-6 then it may be useful to increase the threshold until the warning disappears. Check the “G-matrix eigenvalues” table in the output file for a suitable value.

Branch

GEOMETRY
 Branch {New / Old}
END
Branch
Old / New
Expert option. Specifies which branch of the code to use for making steps. Default the branch of code used depends on the optimization used. Optimization in delocalized coordinates can only be done with the new branch. Optimization in Z-matrix coordinates can only be done with the old branch. In case of Cartesian optimization default the new branch is used, but the old branch can also be used. The new branch can be used (and is default) in geometry optimizations, transition state searches, and in LT. The old branch is default in IRC and NEB, for which one can not use the new branch.

Inithess

GEOMETRY
   Inithess inithessian.file
END
Inithess
Can only be used in combination with the new branch. With this INITHESS subkey it is possible to read a Hessian from a text files. The only argument is the name of the file containing the initial Hessian. The Hessian must be given in full, in non-mass-weighted Cartesian coordinates, and in atomic units (Hartree/Bohr^2).

Hessupd

GEOMETRY
 Hessupd HessUpdate
END
Hessupd
HessUpdate
Specifies how the Hessian matrix is updated, using the gradient values of the current and the previous geometry. The methods available depend on the optimization branch being used. For both the old and new branches, the following options are available: (i) BFGS : Broyden-Fletcher-Goldfarb-Shanno (ii) MS : Murtagh-Sargent (iii) FARKAS : Farkas-Schlegel, Eq. (15) and (16) of Ref. [5] (iv) FARKAS-BOFILL : Farkas-Schlegel-Bofill, Eq. (15) and (14) of Ref. [5] In the old branch, the following extra options are available: (i) DFP : Davidon-Fletcher-Powell (ii) FS : Fletcher switch (iii) HOSHINO : Hoshino In the new branch, the following extra option is available: (i) BAKKEN-HELGAKER: Bakken-Helgaker, see Ref. [6] The default is BFGS for geometry optimizations.

Step

GEOMETRY
 Step {Rad=MaxRadStep} {Angle=MaxAngleStep} {TrustRadius=MaxRadius}
END
Step

Controls that changes in geometry from one cycle to another are not too large:

MaxRadStep
Can only be used in combination with the old branch. An upper bound on changes in Cartesian coordinates or bond lengths, as the case may be. Default: 0.3 angstrom when optimization is carried out in internal coordinates, 0.15 angstrom for Cartesian optimizations.
MaxAngleStep

Can only be used in combination with the old branch. Similarly this option limits changes in bond angles and dihedral angles. Default: 10 degrees. Input for MaxRadStep, MaxAngleStep is in angstrom and degrees respectively, independently of the units used for atomic coordinates input.

Note: Optimization of ring structures carried out in internal (z-matrix) coordinates is sometimes tricky due to the ill-defined last segment of the ring. When problems arise, try Cartesian optimization or consider using smaller limits on the steps (in particular the angles) so as to prevent the program from breaking the ring beyond repair.

MaxRadius
Can only be used in combination with the new branch. By default, the trust radius is set to 0.01 Bohr times the number of atoms with a minimum of 0.2 Bohr. Using the key, the user can override this, setting a constant value. A conservative value is 0.2. A large system (eg 100 atoms) typically needs a larger trust radius (eg 0.8).

DIIS

GEOMETRY
 DIIS {N=NVect} {CYC=Ncyc}
END
DIIS N=NVect CYC=Ncyc
Can only be used in combination with the new branch. NVect is the number of vectors used by the DIIS interpolation method. NCYC is the number of geometry cycles run before the DIIS starts to modify the geometry steps. Default DIIS is used and default N=5 and CYC=0.

Externprogram

GEOMETRY
 Externprogram externprog.exe coord=coords.inp energy=energy.out grad=grads.out
END
Externprogram

Expert option, can only be used in combination with the new branch. Subkey EXTERNPROGRAM has been added to to allow energies and gradients to be calculated by an external program for use in a geometry optimization.

Note that you need to supply information about atomic fragments, such as the basis set, even though these are not actually used in the calculations.

When ADF is ready to perform an energy and gradient calculation, it writes the current cartesian coordinates to the file name given in the input. The format is similar to the ATOMS block in the ADF input file: it has one atom per line, with the element symbol given, followed by the x, y, and z coordinates.

ADF will then run the executable program, and then read in the energy and gradients from the file names given in the input file. The external program is thus responsible for reading the coordinates (in atomic units) written by ADF from file, generating the corresponding energy and gradients (in atomic units), and writing these to the appropriate files. ADF will then take another geometry step, and the process will repeat.

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](1, 2) M. Swart and F.M. Bickelhaupt, Optimization of strong and weak coordinates, International Journal of Quantum Chemistry 106, 2536 (2006)
[5](1, 2) Ö. Farkas and H.B. Schlegel, Methods for optimizing large molecules. Part III. An improved algorithm for geometry optimization using direct inversion in the iterative subspace (GDIIS), Physical Chemistry Chemical Physics 4, 11 (2002)
[6]V. Bakken and T. Helgaker, The efficient optimization of molecular geometries using redundant internal coordinates, Journal of Chemical Physics 117, 9160 (2002)