[hilbert_space] Serious bug in imperative_operator introduced by #686

[krivenko]

Hi everyone,

My collaborators noticed a strange behavior of atom_diag, when energy levels of the same Hamiltonian were different depending on the order of index combinations in the supplied fops list.

The following minimal example demonstrates the problem.

from triqs.operators import c, c_dag, Operator
from triqs.atom_diag import *
from itertools import product

orbs = (1, 2, 3)
spins = ("dn", "up")

# Construct a sum of the pair hopping and spin flip terms from a 3-orbital Kanamori interaction Hamiltonian
H_p = Operator()
H_J = Operator()

for o1, o2 in product(orbs, repeat = 2):
    if o1 == o2: continue
    
    H_p += c_dag(o1, "dn") * c_dag(o1, "up") * c(o2, "up") * c(o2, "dn")
    H_J += c_dag(o1, "dn") * c(o1, "up") * c_dag(o2, "up") * c(o2, "dn")

H = H_p + H_J

fops1 = [(o, s) for (o, s) in product(orbs, spins)]
fops2 = [(o, s) for (s, o) in product(spins, orbs)]

ad1 = AtomDiag(H, fops1)
ad2 = AtomDiag(H, fops2)

for e1, e2 in zip(ad1.energies, ad2.energies):
    print(e1.round(5), e2.round(5))

The output I am getting on my 3.0.x installation is as follows.

[0. 0. 3.] [0. 3. 3.]
[0. 2.] [0. 3. 3.]
[0. 2.] [1. 1. 4.]
[0. 2.] [1. 3.]
[0. 2.] [1. 3.]
[0. 2.] [1. 3.]
[0. 0. 3.] [1. 3.]
[0. 2.] [1. 3.]
[0. 2.] [1. 3.]
[0. 0. 3.] [1. 3.]
[0. 2.] [1. 3.]
[0. 0. 3.] [1. 1. 4.]
[0. 2.] [1. 3.]
[0. 2.] [1. 3.]
[0. 2.] [1. 3.]
[0. 2.] [1. 3.]
[1.] [2.]
[1.] [2.]
[1.] [2.]
...
[1.] [2.]

The energy levels are clearly different (the E = 4.0 level is missing from ad1's output).

I have traced the problem down to imperative_operator's constructor and to PR #686, which broke an assumption that constructor relied upon. More details are provided in the proposed solution #820.

All versions and branches since 2.2 seem to be affected.

[krivenko]

I'd like to stress that this issue is pretty bad. It has the potential to invalidate all CTHYB results for non density-density multiorbital models obtained since the 2.2 release.

[Wentzell]

Dear @krivenko,

Thank you for pointing this out!
This is indeed a serious issue, and I am quite struck how this was able to slip through the tests
of both atom_diag and cthyb.

I will review your proposed solution and will of course back port the fix to the 2.2.x branch.

[Wentzell]

This issue was resolved on all effected branches, c.f. commits 35b968b, babab7d and fee87d7.

[Wentzell]

I have also put in place a simple mechanism (c.f. a54822a) to detect if a given model was effected by the issue.
By setting the environment variable CHECK_ISSUE819 to an arbitrary value with e.g. export CHECK_ISSUE819=1
before starting the calculation, TRIQS will now throw an exception if the model was effected by the issue.