Diagonalizations are important part of (relativistic) quantum chemistry methods.

In the DIRAC suite, there is the diagonalization of the Lowdin matrix, then of the Jz MO matrix in the case of linear symmetry and, finally, of the Fock MO matrix.

In DIRAC, there are few routines to accomplish these tasks:

- non-deafult lapack’ routine DSYEVR (for real symmetric matrixes only)
- eispack’s routine RS (for real symmetric matrixes only)
- own DIRAC routine JACOBI (for real symmetric matrices only)
- own DIRAC routine QJACOBI (for quaternion Hermitian matrices)
- and own DIRAC routine based on Householder diagonalization (for quaternion Hermitian matrices)

For certain cases one has to carefully choose the diagonalization method, depending either on the quality (accuracy) of obtained results - eigenvalues and eigenvectors - or on the time performance.

Note that DSYEVR is not the default diagonalization routine in DIRAC, but can be activated by *./setup -D LAPACK_QDIAG=1 ...*.

The lapack’ routine for diagonalization, DSYEVR, which is fast enough, is giving - for certain input matrices - numerically unstable eigenvectors. We found it in the LINSYM routine when the linear symmetry is turned on.

The real test system, where numerical discrepancies with active DSYEVR are observed, is the \(O_2\) molecule of the *fscc_highspin* test:

```
pam --noarch --inp=hf.inp --mol=O2.mol
```

(Input files for download are from the DIRAC test directory, namely `hf.inp` and
`O2.mol` )

The proper Total SCF energy of this system is,

```
-149.68666145184687
```

while the DSYEVR routine employed in the LINSYM procedure gives different value by about \(10^{-7}\), like (obtained by Miro)

```
-149.68666198782267
```

together with some positron states intruding.

Molecular symmetries lower than linear symmetry (D2h, C2v ...) give proper and identical SCF energies.

By tracing the error to the smallest system where error is clearly showing up we came to the artificial *X2* one-electron linear system
in the smallest 1s1p1d basis:

```
pam --noarch --inp=dc.1el.2fs.inp --mol=X2.1s1p1d.asd.mol --get "Jz_SS_matrix.fermirp2-2"
pam --noarch --inp=dc.1el.2fs.inp --mol=X2.1s1p1d.D2h.mol
```

(Input files for download are `dc.1el.2fs.inp`,
`X2.1s1p1d.asd.mol` and
`X2.1s1p1d.D2h.mol`.)

The affected value is the “Inactive energy (Output from ONEERG)” in the output:

```
-15622.826231066265 (wrong value, linear symmetry)
-15622.826368346343 (correct value, linear symmetry)
-15622.826368346305 (correct value, D2h symmetry)
```

We see differences in energies by about \(10^{-4}\) a.u. what is sufficient for detecting an error.

By tracking down the source of the problem we get to the 9x9 matrix of the *X2.1s1p1d* test system above, *Jz_SS_matrix.fermirp2-2*,
for the diagonalization testing, which is saved in the LINSYM routine.

In the linear symmetry eigenvectors obtained with the lapack’s DSYEVR routine happen to be numerically inaccurate. This can be verified with the abovementioned matrix.

With the standalone DIRAC-dependent diagonalization testing program, *diag.x* (its source code is in the utils directory,
together with the text input file `DIAGONALIZE_TESTS.INP`), we found discrepancies
in numerical accuracy of eigenvectors (U)
obtained with the DSYEVR routine in comparison to other diagonalization methods. Eigenvalues (eps) are unaffected.
A is the original matrix for diagonalization, I is the unit matrix. The norm is defined as sum of absolute values of elements divided by total number of elements.
We split matrix elements into diagonal and off-diagonal.

Values around \(10^{-6} - 10^{-5}\) clearly indicate deviation from the zero threshold for which one would expect numbers cca \(10^{-15} - 10^{-14}\).

By using the testing matrix, *Jz_SS_matrix.fermirp2-2*, we get properly degenerate eigenvalues
with DSYEVR and with other diagonalization methods:

```
1 -1.50000000000001
2 -1.50000000000000
3 -1.50000000000000
4 0.499999999999999
5 0.500000000000000
6 0.500000000000000
7 0.500000000000000
8 0.500000000000002
9 2.50000000000000
```

Now we can proceed to checking the numerical accuracy of obtained eigenvectors.

```
U^{+}*A*U - eps = 0 > norm/diag:0.1684D-14 norm/offdiag:0.6623D-06
U^{+}*U -I = 0 > norm/diag:0.2097D-15 norm/offdiag:0.1325D-05
U*U^{+} -I = 0 > norm/diag:0.6293D-05 norm/offdiag:0.4174D-05
```

```
U^{+}*A*U - eps = 0 > norm/diag:0.1739D-14 norm/offdiag:0.5091D-06
U^{+}*U -I = 0 > norm/diag:0.2097D-15 norm/offdiag:0.1018D-05
U*U^{+} -I = 0 > norm/diag:0.5239D-05 norm/offdiag:0.3407D-05
```

```
U^{+}*A*U - eps = 0 > norm/diag:0.5058D-15 norm/offdiag:0.9620D-16
U^{+}*U -I = 0 > norm/diag:0.6168D-15 norm/offdiag:0.1400D-15
U*U^{+} -I = 0 > norm/diag:0.5674D-15 norm/offdiag:0.1586D-15
```

For public purposes we are offering the standalone and DIRAC independent demonstration program *dsyerv_check.F90*
which does the same checks as the above mentioned (DIRAC dependent) *diag.x* program.
(Source codes for downloading are `dsyerv_check.F90`,
`eispack.F`, and the testing matrix is
`Jz_SS_matrix.fermirp2-2`.)

Results are as follows:

```
U^{+}*A*U - eps ?= 0> norm/diag:0.1783D-14 norm/offdiag:0.1419D-06
U^{+}*U - I ?= 0> norm/diag:0.2591D-15 norm/offdiag:0.2838D-06
U*U^{+} - I ?= 0> norm/diag:0.1047D-05 norm/offdiag:0.7950D-06
```

```
U^{+}*A*U - eps ?= 0> norm/diag:0.1684D-14 norm/offdiag:0.2746D-15
U^{+}*U - I ?= 0> norm/diag:0.2097D-15 norm/offdiag:0.6923D-16
U*U^{+} - I ?= 0> norm/diag:0.2467D-16 norm/offdiag:0.8681D-16
```

It is interesting to observe that obtained Intel+MKL results are in fact correct and
different that those from the *diag.x* DIRAC based standalone program.
This can be attributed to more complicated compiling and linking flags
of the *diag.x* code in the complex DIRAC’s cmake buildup apparatus.

```
U^{+}*A*U - eps ?= 0> norm/diag:0.1665D-14 norm/offdiag:0.9524D-07
U^{+}*U - I ?= 0> norm/diag:0.2097D-15 norm/offdiag:0.1905D-06
U*U^{+} - I ?= 0> norm/diag:0.7649D-06 norm/offdiag:0.5035D-06
```

- The “newest” (in August 2014) lapack 3.5.0 from netlib still suffers with non-orthogonal eigenvectors
- (see this bug report ).

The simplest solution is to keep non-DSYEVR diagonalization routine in the LINSYM part to avoid possibility of highly degenerate values with numerically unstable eigenvectors. The eispack’s RS routine, own DIRAC’s Householder or Jacobi methods do fit for this purpose because they can properly handle eigenvectos of degenerate eigenvalues. The drawback is the time efficiency as the diagonalized matrix is of the Fock matrix size.

More elaborate solution for improving the diagonalization methods for degenerate eigenvalues is pointed out by Luuk Visscher (private communication).