DETCI: incorrect roots and diagonalizer inconsistencies

[bgpeyton]

I recently found an issue when using DETCI to compute roots of LiCN with CIS/6-31G*, which may have uncovered some deeper problems.

the main issue: the lowest-energy solution as predicted by GAMESS (and my own code, which diagonalizes the Hamiltonian by simply calling np.linalg.eigh) is completely skipped by Psi4, despite the following roots being the same. GAMESS predicts the following first six excited state energies (the RHF energy matching to 6 decimals):

GAMESS LiCN/6-31G*
CONVERGED STATE    1 ENERGY=      -99.5348374214
CONVERGED STATE    2 ENERGY=      -99.5315000308
CONVERGED STATE    3 ENERGY=      -99.5315000308
CONVERGED STATE    4 ENERGY=      -99.5258078193
CONVERGED STATE    5 ENERGY=      -99.5258078193
CONVERGED STATE    6 ENERGY=      -99.5095744807

while Psi4 predicts:

Psi4 LiCN/6-31G* (symmetry c1, default [Davidson] solver)
CI ROOT 1 TOTAL ENERGY:-99.5315015582624
CI ROOT 2 TOTAL ENERGY:-99.53150155297872
CI ROOT 3 TOTAL ENERGY:-99.52581092386518
CI ROOT 4 TOTAL ENERGY:-99.52581092379563
CI ROOT 5 TOTAL ENERGY:-99.5096013650753
CI ROOT 6 TOTAL ENERGY:-99.48378441050741

this result doesn't change with or without frozen core, but the results are fine with STO-3G, as well as water with 6-31G*. the results also match GAMESS if I perturb the linear geometry of LiCN, as suggested by @psi-rking , which suggests maybe an issue with finding roots of whatever symmetry the first excited state is; however, adding more roots (or running with/without symmetry) doesn't resolve the issue. in fact, running in c2v produces different energies, but they don't match GAMESS or the symmetry c1 case:

Psi4 LiCN/6-31G* (symmetry c2v, default [Davidson] solver)
CI ROOT 1 TOTAL ENERGY:-99.52581092382692
CI ROOT 2 TOTAL ENERGY:-99.50960136455492
CI ROOT 3 TOTAL ENERGY:-99.47046293582866
CI ROOT 4 TOTAL ENERGY:-99.43820506440449
CI ROOT 5 TOTAL ENERGY:-99.39012366739456
CI ROOT 6 TOTAL ENERGY:-99.37979350854745

maybe it's a preconditioner issue? the Davidson seems to converge in just one step for c1 and c2v. I suppose it's possible the guess space just isn't grabbing roots for all symmetries, but increasing the number of roots to 30 doesn't change the results.

perhaps more disturbing, though, is that the in-core diagonization routine RSP gives yet another set of results:

Psi4 LiCN/6-31G* (symmetry c1, RSP solver)
CI ROOT 1 TOTAL ENERGY:-99.58887523722994
CI ROOT 2 TOTAL ENERGY:-99.56069933848028
CI ROOT 3 TOTAL ENERGY:-99.56069933848008
CI ROOT 4 TOTAL ENERGY:-99.53483755106353
CI ROOT 5 TOTAL ENERGY:-99.53483755106343
CI ROOT 6 TOTAL ENERGY:-99.53418772103414

here, we find 3 roots lower in energy than either GAMESS or the default Psi4 routines, and we also recover the -99.534837 E_h root, but this time as a doubly-degenerate state. the other routines mentioned in the manual
are deprecated kwargs (see also: #2640 ), so I can't test them.

given np.linalg.eigh just uses LAPACK's syevd/heevd, and RSP is just supposed to be LAPACK's dsyev, this result is very perturbing.

I've attached output files for the Davidson and RSP schemes, as well as GAMESS.