RFO and updates

As we said in section 1.3 Rational Function Optimization is a second order method[130,131,245]. It is basically a Newton-Raphson procedure, a Hessian eigenmode following algorithm that avoids the inversion of the Hessian matrix and provides an implicit determination of step length.

Although we already described the mathematics of RFO method in the introduction section 1.3.4.2, now we will give a handful explanation of the method that will stand out the advantages and the difficulties of its practical usage.

RFO is an iterative process that can be outlined as follows:

We first build the Augmented Hessian(AH), which is the Hessian matrix (${\bf B}_k$) with an additional row and an additional column that contains the gradient.

$$AH=\left(\begin{array}{cc} 0 & {\bf g}^T_k \\ {\bf g}_k & {\bf B}_k \end{array}\right)$$ (4.1)

At every iteration $k$ the (N+1) eigenvalue system must be solved

$$\left(\begin{array}{cc} 0 & {\bf g}^T_k \\ {\bf g}_k & {\bf B}_k... ...mbda^{(k)}_\theta {\bf v}^{(k)}_\theta \qquad \forall \quad \theta=1,\ldots,N+1$$ (4.2)

The displacement vector $\Delta {\bf q}_k$, is obtained dividing the first ($\theta=1$) or second ($\theta=2$) eigenvector (depending if we look for a minimum or a saddle point respectively) by its first component

$$\Delta {\bf q}_k=\frac{1}{ v^{(k)}_{1,\theta}}{\bf v}^{'(k)}_\theta$$ (4.3)

The quadratic variation energy $Q(\Delta {\bf q}_k)$ is evaluated as

$$Q(\Delta {\bf q}_k)=\frac{1}{2}\frac{\lambda^{(k)}_\theta}{v^{(k)}_{1,\theta}}$$ (4.4)

It can be shown that as the optimization process converges $v^{(k)}_{1,\theta}$ tends to 1 and $\lambda^{(k)}_\theta$ to zero. This criteria is an additional test of convergence aside from the gradient norm and the change in energy (section 1.3.1).

The $(N+1)$ eigenvalue system (equation 3.2) may be solved completely if we want to check at every step the correct concavity of the PES. If this is not the case, a partial diagonalization that only gives the two lowest eigenpairs will be enough to calculate the displacement (equation 3.3) and proceed with the search^4.2.

Updating formula for the Hessian matrix:
When we are trying to locate a minimum with RFO, the Hessian matrix is updated using a modified form of the BFGS formula[12]. The most general rank-two update Hessian matrix formula is

$${\bf B}_{k+1}={\bf B}_0+\sum_{i=0}^k[{\bf j}_i{\bf u}_i^T +{\bf u... ..._i^T -({\bf j}_i^T \Delta {\bf q}_i){\bf u}_i {\bf u}_i^T ] \qquad k=0,1,\ldots$$ (4.5)

Where ${\bf j}_i={\bf D}_i-{\bf A}_i$, ${\bf D}_i={\bf g}_{i+1}-{\bf g}_i$, ${\bf A}_i={\bf B}_i\Delta {\bf q}_i$, ${\bf u}_i={\bf M}_i\Delta {\bf q}_i / (\Delta {\bf q}^T_i{\bf M}_i\Delta {\bf q}_i)$ and ${\bf B_{k+1}}$ is the approximated Hessian matrix. Different election of the ${\bf M}_i$ matrix leads to different update Hessian matrix formula. In particular, for the BFGS update ${\bf M}_i=a_i{\bf B}_{i+1}+b_i{\bf B}_i$ for some selected positive definite scalars $a_i$ and $b_i$[136].

The proposed modified form of the BFGS expression only differs with respect to the normal BFGS in the calculation of the two scalars $a_i$ and $b_i$. In this modified form these two scalars are evaluated as

$$a_i=\frac{[{\bf A}^T_i {\bf D}_i]^2}{({\bf A}^T_i{\bf A}_i)({\bf D}^T_i){\bf D}_i}$$ (4.6)

$$b_i= 1-a_i$$ (4.7)

Note that both $a_i$ and $b_i$ are positive quantities. The resulting ${\bf B}_{k+1}$ updated Hessian will be positive-definite if both the ${\bf B}_{k}$ matrix and the scalar $\Delta {\bf q}^T_k{\bf D}_k$ are positive-definite as occurs in the normal BFGS formula.

For location of first-order saddle points the Powell formula is used, where in this case the matrix ${\bf M}_i$ contained in equation 3.5 is equal to unit matrix ${\bf I}$.

Acceleration of optimization process: DIIS
We implemented the Direct Inversion of Iterative Space to the gradient vector (GDIIS) as a process to accelerate the convergence in the vicinity of the stationary point. As described in section 1.3.4.3 a combination of the previous gradient vectors is used so as to minimize the length of the current gradient vector.

The GDIIS procedure has only been used in a development stage of the source code. We saw that the quality of a good initial Hessian could improve in a more efficient way the convergence process. So it is not used in the test systems nor in the forthcoming sections of application to enzymatic systems.