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:

1. Murtagh/Sargent (MS):

   $$F^k = F^{k-1} + \frac{ Z^{k-1} (Z^{k-1})^\dagger } { (Z^{k-1})^\dagger dG^{k-1} }$$

   
   

2. Broyden/Fletcher/Goldfarb/Shanno (BFGS):

   $$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}$$

   
   

3. Davidon/Fletcher/Powell (DFP):

   $$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 }$$

   
   

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

$$H^k_{ii} = \sqrt{H^{k-1}_{ii} (h_i + d_i) }$$

with the new estimate $h$ for the diagonal elements obtained by

$$h_i= \frac{\sum_k dG^k_i dq^k_i}{\sum_k (dq^k_i)^2}$$

and the error $d$ obtained by the regression

$$d_i= \frac{\sqrt{ \frac{ \sum_k (dg^k_i)^2} { \sum_k (dq^k_i)^2} - h_i^2 }}{k-2} \:.$$

Another alternative is to use DIIS-like methods: structure optimization by direct inversion in the iterative subspace. (See ref. [29] 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.

$\vert\vert G\vert\vert < \mathrm{thr}$:

BFGS-type update of the hessian and quasi-Newton update of (generalized) coordinates.

References for the algorithms mentioned above: [30,26,31,29,32,33]