5.3 Program Relax

5.3.1 Purpose

relax drives and controls a non-linear optimization procedure to locate the minimum (or a stationary point) of a function $f(x)$. In TURBOMOLE $f$ is always the electronic energy, and the coordinates $x$ will be referred to as general coordinates. They include

cartesian atomic coordinates
internal atomic coordinates
exponents, contraction coefficients and scaling factors of basis functions
a global scaling factor (a common scaling factor for all basis set exponents)

The optimization employs an iterative procedure based on gradients $\nabla f$ of the current and, if available, previous iterations. Various procedures can be applied: steepest descent, Pulay’s DIIS, quasi–Newton, conjugate gradients, as well as combinations of them. relax carries out:

update of general coordinates
update of approximate Hessians if needed
conversion of coordinates (internal $\longleftrightarrow$ Cartesian)

The mode of operation is chosen by the keywords $optimize and $interconversion and the corresponding options, which will be described in the following sections.

5.3.2 Optimization of General Coordinates

After gradients $G^k$ have been calculated for coordinates $q^k$ in optimization cycle $k$, new coordinates (or basis set exponents) $q^{k+1}$ can be obtained from the quasi–Newton update: \[\begin{equation*} q^{k+1} = q^k - F^k G^k \end{equation*}\] where $F^k$ is the inverse of an approximate force constant matrix $H^k$. This method would immediately converge to the equilibrium geometry if $F^k$ would be the inverse of the exact force constant matrix and the force field would be quadratic. In real applications usually none of these requirements is fulfilled. Often only a crude approximation to the force constant matrix $H^k$ is known. Sometimes a unit matrix is employed (which means coordinate update along the negative gradient with all coordinates treated on an equal footing).

The optimization of nuclear coordinates in the space of internal coordinates is the default task performed by relax and does not need to be enabled. Any other optimization task requires explicit specifications in data group $optimize, which takes several possible options:

$optimize options

internal on/off: Structure optimization in internal coordinates.
redundant on/off: Structure optimization in redundant coordinates.
cartesian on/off: Structure optimization in Cartesian coordinates.
basis on/off: Optimization of basis set exponents, contraction coefficients, scaling factors.
global on/off: Optimization of global scaling factor for all basis set exponents.

Note: All options except internal are switched off by default, unless they have been activated explicitly by specifying on.

Some of the options may be used simultaneously, e.g.

internal, basis
internal, global
cartesian, basis

Other options have to be used exclusively, e.g.

internal, cartesian
basis, global

The update of the coordinates may be controlled by special options provided in data group $coordinateupdate which takes as options:

dqmax=real: Maximum total coordinate change (default: 0.3).
interpolate on/off: Calculate coordinate update by inter/extrapolation using coordinates and gradients of the last two optimization cycles (default: interpolate on) if possible.
statistics integer/off: Display optimization statistics for the integer previous optimization cycles. Without integer all available information will be displayed. off suppresses optimization statistics.

The following data blocks are used by program relax:

Input data from gradient programs grad, rdgrad, egrad, rimp2, mpgrad, etc.:

$grad

Cartesian atomic coordinates and their gradients.

$egrad

exponents and scale factors and their gradients.

$globgrad

global scale factor and its gradient.
Input data from force constant program aoforce:

$grad

Cartesian atomic coordinates and their gradients.

$globgrad

global scale factor and its gradient.

$hessian

the force constant matrix in the space of Cartesian coordinates.
Output data from program relax:

$coord

cartesian atomic coordinates.

$basis

exponents and scale factors.

$global

global scale factor.

For structure optimizations the use of (redundant) internal coordinates is recommended, see Section 4. Normally internal coordinates are not used for input or output by the electronic structure programs (dscf, mpgrad, etc.). Instead the coordinates, gradients, etc. are automatically converted to internal coordinates by relax on input and the updated positions of the nuclei are written in Cartesian coordinates to the data group $coord. Details are explained in the following sections.

5.3.3 Force Constant Update Algorithms

In a Newton-type geometry update procedure often only a crude approximation to the force constant matrix $H^k$ is available. What can be done then is to update $F^k=(H^k)^{-1}$ in each iteration using information about previous coordinates and gradients. This constitutes the quasi–Newton or variable metric methods of which there are a few variants:

Murtagh/Sargent (MS): \[\begin{equation*} F^k = F^{k-1} + \frac{ Z^{k-1} (Z^{k-1})^\dagger } { (Z^{k-1})^\dagger dG^{k-1} } \end{equation*}\]
Broyden/Fletcher/Goldfarb/Shanno (BFGS): \[\begin{equation*} F^k = F^{k-1} + \frac{ S (dq^{k-1})^\dagger dq^{k-1} - dq^{k-1} (dG^{k-1})^\dagger F^{k-1} - F^{k-1} dG^{k-1} (dq^{k-1})^\dagger } {S1} \end{equation*}\]
Davidon/Fletcher/Powell (DFP): \[\begin{equation*} F^k = F^{k-1} + \frac{(dq^{k-1})^\dagger dq^{k-1} } {S1} - \frac{F^{k-1} dG^{k-1} (dG^{k-1})^\dagger F^{k-1} } {(S-1) S1 } \end{equation*}\]
combined method (BFGS/DFP): If $S1 < (S-1) S1$ and $S1 > 0$ perform DFP update, otherwise BFGS.

The meaning of the symbols above is as follows:

$F^k = (H^k)^{-1}$: approximate inverse force constant matrix in the k-th iteration.s
$q^k$: general coordinates in the k-th iteration.
$G^k$: gradients in the k-th iteration.
$dq^{k-1} = q^k - q^{k-1}$
$dg^{k-1} = g^k - g^{k-1}$
$Z^{k-1} = dq^{k-1} - F^{k-1} dG^{k-1}$
$S1 = (dq^{k-1})^\dagger dg^{k-1}$
$S = 1 + ((dg^{k-1})^\dagger F^{k-1} dG^{k-1} )/ (S1)$

An alternative is to use update algorithms for the hessian $H^k$ itself:

Ehrig, Ahlrichs : Diagonal update for the hessian by means of a least squares fit \[\begin{equation*} H^k_{ii} = \sqrt{H^{k-1}_{ii} (h_i + d_i) } \end{equation*}\] with the new estimate $h$ for the diagonal elements obtained by \[\begin{equation*} h_i= \frac{\sum_k dG^k_i dq^k_i}{\sum_k (dq^k_i)^2} \end{equation*}\] and the error $d$ obtained by the regression \[\begin{equation*} d_i= \frac{\sqrt{ \frac{ \sum_k (dg^k_i)^2} { \sum_k (dq^k_i)^2} - h_i^2 }}{k-2} \:. \end{equation*}\] Another alternative is to use DIIS-like methods: structure optimization by direct inversion in the iterative subspace. (See ref. [59] for the description of the algorithm). The DIIS procedure can often be applied with good success, using static or updated force constant matrices.

Any of the algorithms mentioned above may be chosen. Recommended is the macro option ahlrichs, which leads to the following actions ($n$ is the maximum number of structures to be included for the update, default is $n=4$):

$\mathrm{ncycles} < n$:: geometry update by inter/extrapolation using the last 2 geometries.
$\mathrm{ncycles} \ge n$:: diagonal update for the hessian as described above; DIIS–like update for the geometry.
$||G|| < \mathrm{thr}$:: BFGS-type update of the hessian and quasi–Newton update of (generalized) coordinates.

References for the algorithms mentioned above: [60, 56, 61, 59, 62, 63]

5.3.4 Definition of Internal Coordinates

If structure optimizations are to be performed in the space of internal coordinates ($optimize internal, is the default setting), appropriate internal coordinate definitions have to be provided on data block $intdef. The types available and their definitions are described in Section 4.1.2. For recommendations about the choice of internal coordinates consult ref. [49]. Nevertheless the structure of $intdef will shortly be described. The syntax is (in free format):

 1   k   1.00000000   bend  1  2  3    val=1.9500  fdiag=.6666

The first items have been explained in Chapter 4.

Two additional items val=real, fdiag=real may be supplied for special purposes:

val=: serves for the input of values for internal coordinates for the interconversion internal $\rightarrow$ cartesian coordinates; it will be read in by relax if the flag for interconversion of coordinates has been activated ($interconversion on ), or by the interactive input program define within the geometry specification menu.
fdiag=: serves for the input of (diagonal) force constants for the individual internal coordinates to initialize $forceapprox.

5.3.5 Structure Optimizations Using Internal Coordinates

This is the default task of relax ($optimize internal on does not need to be specified!) You need as input the data groups :

$grad: Cartesian coordinates and gradients as provided and accumulated in subsequent optimization cycles by the programs grad, or rdgrad etc.
$intdef: definitions of internal coordinates.
$redundant: definitions of redundant coordinates.

Output will be the updated coordinates on $coord and the updated force constant matrix on $forceapprox. If any non-default force constant update option has been chosen, relax increments its counting variables <numgeo>, <numpul> within command keyword $forceupdate. If the approximate force constant has been initialized ($forceinit on ) relax switches the initialization flag to $forceinit off. Refer also to the general documentation of TURBOMOLE. It is recommended to check correctness of your definition of internal coordinates:

Calculate their values for your Cartesian start coordinates using the relax program (see Section 5.3.11) or within a define session.
Have a look at the eigenvectors of the $\mathbf{BmB}^\dagger$-matrix. Set some ‘?’ behind keyword $intdef, if there are any eigenvalues close to zero ($<\!10^{-2}$ is to be considered bad for small molecules, but there is no general rule) check those internal coordinates for consistency which contribute to the corresponding eigenvector(s)!

5.3.6 Structure Optimization in Cartesian Coordinates

For this task you have to specify:

$optimize              
  cartesian on    
  internal  off

These lines switch on the non-default optimization in Cartesian coordinates and switch off the optimization in internal coordinates (this has to be done explicitly!). As input data groups you need only $grad as provided by on of the gradient programs. For the first coordinate update an approximate force constant matrix is needed in data group $forceapprox. Output will be the updated coordinates on $coord, and the updated force constant matrix on $forceapprox.

The coordinates for any single atom can be fixed by placing an ’f’ in the third to eighth column of the chemical symbol/flag group. As an example, the following coordinates specify acetone with a fixed carbonyl group:

$coord
    2.02693271108611      2.03672551266230      0.00000000000000      c
    1.08247228252865     -0.68857387733323      0.00000000000000      c f
    2.53154870318830     -2.48171472134488      0.00000000000000      o     f
   -1.78063790034738     -1.04586399389434      0.00000000000000      c
   -2.64348282517094     -0.13141435997713      1.68855816889786      h
   -2.23779643042546     -3.09026673535431      0.00000000000000      h
   -2.64348282517094     -0.13141435997713     -1.68855816889786      h
    1.31008893646566      3.07002878668872      1.68840815751978      h
    1.31008893646566      3.07002878668872     -1.68840815751978      h
    4.12184425921830      2.06288409251899      0.00000000000000      h
$end

5.3.7 Optimization of Basis Sets

For this task you have to specify:

$optimize              
  basis     on   
  internal  off

This example would perform only a basis set optimization without accompanying geometry optimization. It is possible, of course, to optimize both simultaneously: Just leave out the last line of the example (internal off). Input data groups are:

$egrad: Basis set exponents, contraction coefficients, scaling factors and their respective gradients as provided and accumulated in subsequent optimization cycles by one of the programs grad or mpgrad, if $drvopt basis on has been set.
$basis: Description of basis sets used, see Section 4.2.

Output will be the updated basis on $basis, and the updated force constant matrix on $forceapprox.

For an example, see Section 26.5.

Note that the optimization of exponents and scaling factors is restricted to the Hartree–Fock level of theory. Optimization of contraction coefficients can be carried with Hartree–Fock or DFT.

5.3.8 Simultaneous Optimization of Basis Set and Structure

The optimization of geometry and basis set may be performed simultaneously and requires the specification of:

$optimize              
  internal  on   (or: cartesian on)
  basis     on

and needs as input data groups $grad and $egrad. Output will be on $coord, $basis, also on $forceapprox (updated).

5.3.9 Optimization of Structure and a Global Scaling Factor

Optimization of a global scaling factor is usually not performed in geometry optimizations. It is a special feature for special applications by even more special users. As reference see [64].

To optimize the structure and a global scaling factor specify:

$optimize              
  internal  on   (or: cartesian on)
  global    on

You need as input data groups $grad and $globgrad, the latter contains the global scaling factors and their gradients accumulated in all optimization cycles. Output will be on $coord, $global, also on $forceapprox (updated). Note that for optimization of a global scaling factor a larger initial force constant element is recommended (about $10.0$).

5.3.10 Conversion from Internal to Cartesian Coordinates

Due to translational and rotational degrees of freedom and the non-linear dependence of internal coordinates upon Cartesian coordinates, there is no unique set of Cartesian coordinates for a given set of internal coordinates. Therefore an iterative procedure is employed to calculate the next local solution for a given Cartesian start coordinates. This task may be performed using the relax program, but it is much easier done within a define session.

5.3.11 Conversion of Cartesian Coordinates, Gradients and Force Constants to Internals

To perform this tasks, you have to activate the interconversion mode by

$interconversion on
   cartesian --> internal   coordinate gradient hessian

Note that any combination of the three options showed is allowed! The default value is coordinate, the two other have to be switched on explicitly if desired.

You need as input data groups:

intdef: Definitions of (redundant) internal coordinates
coord: Cartesian coordinates (for option ‘coordinate’)
grad: Cartesian coordinates and gradients as provided and accumulated in subsequent optimization cycles by the various gradient programs (for coordinate and gradient)
hessian: Analytical force constant matrix (as provided by the force constant program aoforce) (only if option hessian is specified). The data group $hessian (projected) may be used alternatively for this purpose.

All output will be written to the screen except for option hessian (output to data group $forceapprox)

5.3.12 The m-Matrix

The m-matrix serves to fix position and orientation of your molecule during geometry optimizations. It cannot be used to fix internal coordinates! The m-matrix is a diagonal matrix of dimension $3n^2$ (where $n$ is the number of atoms). Normally m will be initialized as a unit matrix by relax. As an example consider you want to restrict an atom to the xy-plane. You then set the m(z)–matrix element for this atom to zero. You can use at most six zero m-matrix diagonals (for linear molecules only five)—corresponding to translational and rotational degrees of freedom. Note that the condition of the $\mathbf{BmB}^\dagger$-matrix can get worse if positional restrictions are applied to the m-matrix. m-matrix elements violating the molecular point group symmetry will be reset to one. Non-default settings for m-matrix diagonals of selected atoms have to be specified within data group $m-matrix as:

$m-matrix
     1   0.0  0.0  0.0
    10   1.0  0.0  0.0
    11   1.0  1.0  0.0

5.3.13 Initialization of Force Constant Matrices

The most simple initial hessian is a unit matrix. However, better choices are preferable. For structure optimizations using internal coordinates you may use structural information to set up a diagonal force constant matrix with elements chosen in accord to the softness or stiffness of the individual modes. For detailed information refer to ref. [62]. For optimization of basis set parameters less information is available. When neither data block $forceapprox is available nor $forceinit on is set, the force constant matrix will be initialized as a unit matrix. Specifying the force constant initialization key $forceinit on diag=... will lead to:

diag=real

Initialization with real as diagonal elements.

diag=default

Initial force constant diagonals will be assigned the following default values:

internal coordinates	:	stretches	0.50
		angles	0.20
scaling factors	:	s,p	1.50
		d	3.00
exponents	:	uncontracted	0.15
		contracted	10.00
contraction coefficients	:		100.00
global scaling factor	:		15.00
cartesian force constants	:		0.50

diag=individual

Initial force constant diagonals will be taken from
$intdef fdiag=... or
$global fdiag=...
Similar initialization modes are NOT supported for geometry optimization in Cartesian space and for the optimization of basis set parameters!

carthess

Data group $hessian (projected) is used.

5.3.14 Look at Results

The energy file includes the total energy of all cycles of a structure optimization completed so far. To get a display of energies and gradients use the UNIX command grep cycle gradient which yields, e. g. H$_2$O.

cycle =      1    SCF energy =    -76.3432480651   |dE/dxyz| =  0.124274
cycle =      2    SCF energy =    -76.3575482860   |dE/dxyz| =  0.082663
cycle =      3    SCF energy =    -76.3626983371   |dE/dxyz| =  0.033998
cycle =      4    SCF energy =    -76.3633251080   |dE/dxyz| =  0.016404
cycle =      5    SCF energy =    -76.3634291559   |dE/dxyz| =  0.010640
cycle =      6    SCF energy =    -76.3634910117   |dE/dxyz| =  0.000730

This should be self-evident. To see the current—or, if the optimization is converged, the final—atomic distances use the tool dist. Bond angles, torsional angles etc. are obtained with the tools bend, tors, outp, etc. In the file gradient are the collected Cartesian coordinates and corresponding gradients of all cycles. The values of the general coordinates and corresponding gradients are an output of relax written to job.<cycle> of job.last within jobex. To look at this search for ‘Optimization statistics’ in job.last or job.<cycle>.