Thanks for the detailed explanation, which I believe is very helpful to whoever tries to construct IQMPO by hand. To make sure I do understand what's going on, take the Heisenberg.h as an example (again), where the virtual indices QN is defined as
for(int l = 0; l <= N_; ++l)
q0.at(l) = Index(nameint("q0_",l),3);
qP.at(l) = Index(nameint("qP_",l),1);
qM.at(l) = Index(nameint("qM_",l),1);
links.at(l) = IQIndex(nameint("hl",l),
Then the QN at (row(2), col(4)) is 0+(-2)=-2.
On the other hand, one term in the Hamiltonian is
W += sites_.op("Sp",n) * row(2) * col(4) * J_/2;
where "Sp" has QN +2 as it raises the spin by 1. And -2 and 2 add to be 0. Is this the correct understanding? Thanks!