Wrong MPO when using removeQN for Electron SiteSet

[mingruyang]

I want to do a time evolution which still conserve the quantum numbers but also measure a quantity <(Sy_tot)^2> which violate the Sz conservation. To measure <(Sy_tot)^2> for an initial state with all electrons (half-filling) polarized to the +x direction, I did the following:

int N = 36;

auto sites = Electron(N);

auto is = inds(sites);
for(int i = 1; i <= length(is); ++i)
    is(i) = removeQNs(is(i));
auto sitesQN = SpinHalf(is);

auto ampo = AutoMPO(sitesQN);
for(int i = 1; i <= N; ++i)
    {
    for(int j = 1; j <= N; ++j)
        {
        ampo +=  -0.25,"S+",i,"S+",j;
        ampo +=   0.25,"S+",i,"S-",j;
        ampo +=   0.25,"S-",i,"S+",j;
        ampo +=  -0.25,"S-",i,"S-",j;
        }
    }
auto SySy = toMPO(ampo);

auto state = InitState(sites);
for(int i = 1; i <= N; ++i)
    {
    state.set(i,"Emp");
    }
auto psi1 = MPS(state);

for(int i = N; i >= 1; --i)
    {
    for(int j = 1; j <= i-1; ++j)
        {
        psi1.position(j);
        psi1.Aref(j) = noPrime(op(sites,"F",j)*psi1.A(j));
        }
    psi1.position(i);

    auto ind = removeQNs(sites(i));
    auto indP = prime(ind);
    auto R1 = ITensor(ind, indP);
    R1.set(ind(1),indP(2), 1./sqrt(2));
    R1.set(ind(1),indP(3), 1./sqrt(2));
    R1.set(ind(2),indP(4),-1./sqrt(2));
    R1.set(ind(3),indP(4), 1./sqrt(2));
    psi1.Aref(i) = noPrime(R1*psi1.A(i));
    }

auto sysy = real(innerC(psi1,SySy,psi1));
printfln("sysy=%.10f",sysy);

I will get sysy=4.5, which is wrong.

But if I use instead

auto sites = Electron(N,{"ConserveSz",false});

auto ampo = AutoMPO(sites);
for(int i = 1; i <= N; ++i)
    {
    for(int j = 1; j <= N; ++j)
        {
        ampo +=  -0.25,"S+",i,"S+",j;
        ampo +=   0.25,"S+",i,"S-",j;
        ampo +=   0.25,"S-",i,"S+",j;
        ampo +=  -0.25,"S-",i,"S-",j;
        }
    }
auto SySy = toMPO(ampo);

I will get sysy=9, which is correct.

As mentioned in #311, I guess the wrong result of the first example is because removeQNs also remove the Fermion parity for the Electron SiteSet. Using it the same way, I didn't get wrong results for SpinHalf SiteSet.

Since I do want a SiteSet which conserve Sz so that the time evolution can be faster, and I just want to measure a quantity which breaks the Sz, I was wondering if there is anyway to only remove the Sz QN in this case.

[emstoudenmire]

Thanks for the bug report! I'll look into it and agree it could be the way that removeQNs works that's causing the problem. It's a little confusing here because I don't see where any of the steps would be using the quantum numbers in a way that would change the entries of the tensors, but I'll just have to step through it and put checks along the way to see. AutoMPO does do some complicated things with fermions that could be the culprit.

[mingruyang]

If I just use

auto sites = Electron(N);

auto ampo = AutoMPO(sites);
for(int i = 1; i <= N; ++i)
    {
    for(int j = 1; j <= N; ++j)
        {
        ampo +=  -0.25,"S+",i,"S+",j;
        ampo +=   0.25,"S+",i,"S-",j;
        ampo +=   0.25,"S-",i,"S+",j;
        ampo +=  -0.25,"S-",i,"S-",j;
        }
    }
auto SySy = toMPO(ampo);

I will get an error message as follows:

I = (dim=4|id=876|"l=1,Link") <Out>
  1: 2 QN()
  2: 1 QN({"Nf",0,-1},{"Sz",-2})
  3: 1 QN({"Nf",0,-1},{"Sz",2})
Q = QN({"Nf",0,-1},{"Sz",-4})
From line 683, file index.cc

Index does not contain given QN block.

Index does not contain given QN block.
Aborted (core dumped)