Lanczos

Theory

Ref: [98]

A different way of diagonalizing the Hessian matrix comes from the Lanczos procedure. It is easy to generalize eq.14 to a collective displacement of atoms,

$$H^{(1)}_{PP} = \sum_{R,\alpha}u^R_\alpha \frac \partial {\partial R_\alpha} \; \left[\sum_i \vert\text{\tt P}^R_i\rangle \langle\text{\tt P}^R_i\vert \right]$$ (19)

$$= \sum_{R,\alpha} u^R_\alpha\sum_i\;\left[ \left[\frac\partial{\par... ...[\frac\partial{\partial R_\alpha}\; \langle\text{\tt P}^R_i\vert\right] \right]$$ (20)

In this way the information contained in the Hessian matrix, A, can be used to compute terms of the type

$$\ensuremath{\mathbf{A}}\cdot\ensuremath{\mathbf{w}}= \left( \begi... ...}} {\partial{{R_{Nz}}} \partial{{\ensuremath{\mathbf{w}}}}} }\end{array}\right)$$ (21)

where w is the collective displacement. This is the building block for Lanczos diagonalization, that is performed by iterative application of the symmetric matrix to be diagonalized, over a set of vectors, known as Lanczos vectors, according to the scheme

 $&bull#bullet;$

$\ensuremath{\mathbf{r}}_0=\ensuremath{\mathbf{q}}_1; \beta_0=1; \ensuremath{\mathbf{q}}_0=0; k=0$

WHILE

$\beta_k\neq 0$

 $&bull#bullet;$

$k=k+1$

 $&bull#bullet;$

$\alpha_k=\ensuremath{\mathbf{q}}_k^T\ensuremath{\mathbf{A}}\ensuremath{\mathbf{q}}_k$

 $&bull#bullet;$

$\ensuremath{\mathbf{r}}_k=\ensuremath{\mathbf{A}}\ensuremath{\mathbf{q}}_k-\alpha_k\ensuremath{\mathbf{q}}_k-\beta_{k-1}\ensuremath{\mathbf{q}}_{k-1}$

 $&bull#bullet;$

$\beta_k=\parallel\ensuremath{\mathbf{r}}_k\parallel_2$

 $&bull#bullet;$

$\ensuremath{\mathbf{q}}_{k+1}=\ensuremath{\mathbf{r}}_k/\beta_k$

 $&bull#bullet;$

Orthogonalization, $\ensuremath{\mathbf{q}}_{k+1}\perp\left\{\ensuremath{\mathbf{q}}_1,\cdots,\ensuremath{\mathbf{q}}_k\right\}$

 $&bull#bullet;$

Diagonalization of $\ensuremath{\mathbf{T}}_k$

END WHILE

The matrix $\mathbf{A}$ is thus projected onto a k$\times$k subspace, in a tridiagonal form, $\ensuremath{\mathbf{T}}_k$. $\mathbf{T}$'s eigenvalues are approximations to $\mathbf{A}$'s, while the eigenvectors are brought in the n$\times$n space by means of the orthogonal matrix $\ensuremath{\mathbf{Q}}=[\ensuremath{\mathbf{q}}_1,\ensuremath{\mathbf{q}}_2,\dots,\ensuremath{\mathbf{q}}_k]$. The advantage with respect to the usual diagonalization schemes is that the highest and the lowest end of the spectrum tend to converge before the ionic degrees of freedom are fully explored.

The projection of the trivial eigenmodes is performed by eliminating their contribution from every Lanczos vector at the orthogonalization step reported in the algorithm above,

$$\ensuremath{\mathbf{q}}_k=\ensuremath{\mathbf{q}}_k-\sum_j\left((... ...{\mathbf{q}}_k^T\cdot\ensuremath{\mathbf{s}}_j)\ensuremath{\mathbf{s}}_j\right)$$ (22)

This procedure seems slightly different from that of the PHONON case, but they are perfectly equivalent. It is controlled by the same keyword DISCARD with the same arguments.

Lanczos input

The input for the Lanczos routine is given by the keyword LANCZOS plus some optional arguments and a mandatory following line with numerical data. Plese note that this keyword is analogous to that for the Lanczos diagonalization of the electronic degrees of freedom. Given its presence in the &RESP-section there should be no ambiguity.

• LANCZOS [CONTINUE,DETAILS]; the keyword LANCZOS simply activates the diagonalization procedure. At every cycle the program writes a file, LANCZOS_CONTINUE.1 which contains the information about the dimension of the calculation, the number of iteration, the elements of the Lanczos vectors and the diagonal and subdiagonal elements of the $\mathbf{T}$ matrix.
  • With the option CONTINUE one restart a previous job, asking the program to read the needed information from the file LANCZOS_CONTINUE; if the program does not find this file, it stops. This option is to be used when the calculation that one wants to perform does not fit in the time of a queue slot.
  • The argument DETAILS prints a lot of information about the procedure at the end of the output. It is intended for debugging purposes.
  The subsequent line is mandatory, and contains three numerical parameters;
  • Lanczos_dimension; it is the dimension of the vibrational degrees of freedom to be explored. Normally it is 3$\times$N$_{\rm at}$, -3 if you are eliminating the translations or -6 if you are eliminating also the rotations(vide infra). It is possible to use lower values. Higher are non sense.
  • no. of iterations; it is the number of cycles that are to be performed during the actual calculation. It must be $\leq$ to Lanczos_dimension; the program checks and stops if it is higher. In case of a continued job, the program checks if the given number of iterations + the index of the iteration from which it restarts is within the limit of Lanczos_dimension; if not it resets the number to fit the dimension, printing a warning on std. output.
  • conv_threshold; it is the threshold under which an eigenvector is to be considered as converged. It is the component of the eigenvector over the last Lanczos vector. The lower it is, the lower is the information about this mode which has still to be explored.
• DISCARD {PARTIAL,TOTAL,OFF,LINEAR}; this keyword controls the elimination of the trivial eigenmodes (i.e. translations and rotations) from the eigenvector calculations. It works both for PHONON and LANCZOS. Omitting it is equivalent to DISCARD PARTIAL. When using it, the choice of one of the arguments is mandatory; the program stops otherwise. They are;
  • PARTIAL; only the translational degrees of freedom are eliminated (useful for crystals). This is the default.
  • TOTAL; both translational and rotational degrees of freedom are eliminated.
  • OFF; the option is disactivated and no projection is performed.

Lanczos output In the output, all the informations about the perturbed wavefunction optimization are reported, just like in the PHONON case. The differences are in the form of the perturbation and in the eigenmodes information coming from the diagonalization of $\ensuremath{\mathbf{T}}_k$.

The report of the perturbation for a water dimer reads like;

 **********************   perturbations    **********************

 **** atom=    1     O    .04171821      .07086763      .07833475
 **** atom=    2     O    .03615521      .08499767      .07082773
 **** atom=    3     H    .29222111      .13357862      .09032954
 **** atom=    4     H    .33764012      .15750912      .25491923
 **** atom=    5     H    .06426954      .26020430      .01822161
 **** atom=    6     H    .15765937      .27370013      .29183999
 cpu time for wavefunction initialization:          89.22 seconds

where at every atom corresponds the x,y,z displacemnts applied to calculate the perturbed wavefunctions.

At every iteration information about the elements of the $\ensuremath{\mathbf{T}}_k$ matrix, alpha's and beta's, are reported. Here we are at the end of the calculation. Note the numerical zero in the final +1 beta value. Then the spectrum in a.u. is reported, together with the convergence information. At the last iteration there are vectors which are not ``converged''. But this comes only from the definition, since some of the eigenvectors must have a component over the last Lanczos vectors. Following there are the familiar eigenvalues in cm$^{-1}$.

 ****************************************************************
  *+* L2 norm[n[i+1]]       =  0.100262974102936016E-02
  *=*   overlap: alpha[ 12 ] =  0.982356847531824437E-03
  *=*   off-diag: beta[ 13 ] =  0.300011777832112837E-19
  *=*   norm:               =  0.300011777832112837E-19
 ****************************************************************
 *** SPECTRUM, run 12 :
 *** eigenvalue  12 =     .4855525       (converged: .000000).
 *** eigenvalue  11 =     .4785309       (converged: .000000).
 *** eigenvalue  10 =     .4605248       (converged: .000000).
 *** eigenvalue   9 =     .4303958       (converged: .000000).
 *** eigenvalue   8 =     .0973733       (converged: .000000).
 *** eigenvalue   7 =     .0943757       (converged: .000000).
 *** eigenvalue   6 =     .0141633       (converged: .000004).
 *** eigenvalue   5 =     .0046807   (NOT converged: .004079).
 *** eigenvalue   4 =     .0020023   (NOT converged: .010100).
 *** eigenvalue   3 =     .0010310   (NOT converged: .313417).
 *** eigenvalue   2 =     .0009886   (NOT converged: .933851).
 *** eigenvalue   1 =     .0006303   (NOT converged: .171971).

 ****************************************************************
  harmonic frequencies [cm**-1]:

        129.0591        161.6295        165.0552        230.0241
        351.6915        611.7668       1579.1898       1604.0732
       3372.3938       3488.4362       3555.9794       3581.9733

 ****************************************************************