In this section we recover the RFO equations of section 3.1 and design a method to be used for thousands of atoms. Three main problems arise when working with a big Hessian, that is, the initial Hessian calculation, its diagonalization at every step of the process and its storage. The initial Hessian computational effort is related to the method at which the energy is calculated while the two other issues are only related to the size of the matrix.
The reason why standard Newton methods are not used in systems with a great number of atoms is that computational difficulties are encountered in a big Hessian manipulation. In this chapter we have seen how micro-iterative methods overcome this problem but we must keep in mind that an expensive QM method or an increasing of the core size will provoke these computational problems to reappear.
As we said in the introductory section 1.3.7, memory problems arise in matrices with dimension of 10 000. An example that shows how the computational requirements increase with the size is given in the table below. In this table it is shown the time required to extract the lowest five eigenpairs along with the time needed to calculate the Numerical Hessian at different dimensions of a matrix.
Dimension | Partial diagonalization | Hessian calculation |
99 | 0.01 | 2291 |
285 | 0.16 | 6719 |
420 | 1.1 | 9878 |
978 | 17.3 | 23071 |
3894 | 1195.0 | 92196 |
The given matrix dimension should properly be indicated as | ||
Time given in seconds on a Pentium IV 2.0 GHz |
We do not need bigger Hessian matrices in enzymatic reactions:
One of the conclusions in the previous sections is that a core size higher than a thousand of atoms is not useful.
The storage of a matrix of few hundreds of atoms is not a real problem in nowadays computers.
In addition to the storage, the full diagonalization can be avoided
with a standard partial diagonalization which provides the two lowest eigenpairs needed to calculate
the displacement in the RFO scheme.
Depending on the QM level a compromise must be found between the quality of the initial Hessian and the steps required to converge. That is, a very cheap second derivatives matrix could be used but this would imply a bigger number of steps. Since the number of steps required to converge with a bad initial Hessian can increase very rapidly, our recommendation would be that the higher the quality of the Hessian the better. Of course in a transition state search this last statement is crucial. It has also been shown[151,161] that an adequate initial Hessian eliminates the known problems of coupling between Cartesian coordinates.
We still need a general algorithm:
The above discussion is valid for transition state search of a chemical step in enzymatic reactions.
However, we must take into account that trying to find a general second-order algorithm capable to move a huge
number of atoms is still a very active area.
This ideal algorithm would be needed to minimize a very big number of atoms or to locate saddle points when many
atoms are involved in the transition (e.g. phase transition in solid state[142]).
In this case, the problem of storage of a big Hessian matrix emerges as the most important bottleneck in the algorithm.
Keeping this in mind we have designed a second order method that avoids the storage of the Hessian matrix. It is based on a limited-memory update combined by a Lanczos-type iterative diagonalization that permits to extract the displacement eigenpair at very low memory cost.
This procedure has already been proposed by Bofill and co-workers[155] but in that case only small systems were tested and the usage of internal coordinates may have notable implications that we will comment immediately. Besides, here we analyze the shape of different initial Hessians and we will alert to the problems related to the convergence in the iterative diagonalization.