Adding IQMPOs with different Arrow structure

Question

Hi there,
I am still struggling with adding IQMPOs (IQMPSs) and I hope to get some help here.

A couple of weeks ago I already pointed out that the plusEq function gives segfaults for any IQMPS and IQMPO I tried. After some research I figured out that this is due to the plussers(IQIndex, IQIndex, IQIndex, IQTensor, IQTensor) function. I rewrote it and now I come up with

bool
hasQN(const IQIndex& I, const QN& qn)
    {
    for(auto& iq : I)
         if(iq.qn == qn) return true;
    return false;
    }

void
plussers(IQIndex const& l1, IQIndex const& l2,
         IQIndex& sumind,
         IQTensor& first, IQTensor& second)
    {

    vector<IndexQN> iq;
    vector<IQIndex> map;

    for(IndexQN const& x : l1)
        {
        if(hasQN(l2, x.qn))
            {
            Index y = findByQN(l2, x.qn);
            Index jj(x.index.rawname() + y.rawname(),x.m() + y.m(),x.type());
            iq.push_back(IndexQN(jj, x.qn));
            for (int i = 1; i <= x.m(); i++) { map.push_back(l1); }
            for (int i = 1; i <= y.m(); i++) { map.push_back(l2); }
            }
        else
            {
            Index jj(x.index.rawname(),x.m(),x.type());
            iq.push_back(IndexQN(jj, x.qn));
            for (int i = 1; i <= x.m(); i++) { map.push_back(l1); }
            }
        }
    for(IndexQN const& y : l2)
        {
        if(!hasQN(l1, y.qn))
            {
            Index jj(y.index.rawname(),y.m(),y.type());
            iq.push_back(IndexQN(jj, y.qn));
            for (int i = 1; i <= y.m(); i++) { map.push_back(l2); }
            }
        }

    sumind = IQIndex(sumind.rawname(),std::move(iq),sumind.dir(),sumind.primeLevel());

    assert(sumind.m() == l1.m() + l2.m());
    first = IQTensor(dag(l1), sumind);
    int count = 1;
    for (int i = 1; i <= sumind.m(); i++)
        {
        if (map[i-1] == l1)
            {
            first.set(l1(count), sumind(i), 1.0);
            count++;
            }
        }

    second = IQTensor(dag(l2), sumind);
    count = 1;
    for (int i = 1; i <= sumind.m(); i++)
        {
        if (map[i-1] == l2)
            {
            second.set(l2(count), sumind(i), 1.0);
            count++;
            }
        }
    }

which works quite well as long as all links in the IQMPOs (IQMPSs) have corresponding arrow directions. I mean, when adding IQMPO A and IQMPO B, the first link of A needs to have the same arrow as the first link of B; the second link of A needs to have the same arrow as the second link of B....

But as soon as the n-th Link of A and B differ, the code breaks in addAssumeOrth when adding the parts together. This is due to a QN issue I don't know how to solve.

I hope somebody can give me hint on how to fix this.

I also would be happy about some theory. How are MPOs and MPSs with different arrow structure supossed to be added? I hardy find anything on this in the literature. I tried repositioning of the MPOs (MPSs) but that change the arrow direction in my case (might this be a bug?).

Best,
Lars

Comments

I might got a step in the right direction. If I call A.orthogonalize() and checkQNs(A) (B.orthogonalize() and checkQNs(B)) I get an error, so maybe there is something with the MPOs I'm trying to add wrong. Hope someone can confirm this.

---

Hi Lars,
As you can see I'm finally finding some time to fix more bugs. In the above question there seem to be many different possible error sources. In your comment you say that perhaps the IQMPO's you were adding were ill-formed, in the sense of not have a well defined quantum number? If so I could see how that could break some code later on, such as in plussers and elsewhere.

Can you confirm if it's working now? If not I'd be happy to take a look but I would need to use the same MPO's you were working with i.e. a minimal sample code I could debug from.

Thanks,
Miles

---

Oops, ok I just found the sample code you posted on Github. I will try it out and get back to you with what I find.

---

Hi Miles,

thanks for working on my bug reports. I am not sure whether this is still a bug, my confusion mostly came from the strange behaviour  of the .orthogonalize() function I reported here http://itensor.org/support/241/strange-orthogonalize-behaviour-when-operators-multiplied. And I see you are working on that. So I it would be best to close this.

However, what I think is wrong with the code is on the one hand an infinite recursion in the plusEq function for IQMPOs and something with the plussers function. Would be nice if you can confirm that, but I can also provide examples if you want.

I fixed both errors in my https://github.com/lfrahm/ITensor/tree/cmplmntryoprtrs fork. You might want to have a look at the mpsalgs.cc file. There you find my latest plussers function. It turned out the one I wrote down in my question above had a bug. The one on github works fine so far.

I hope this helps, let me know if I can help you with anything.
Best,
Lars

Answer 1

Hi Lars,
Thanks for your patience. I just pushed a set of patches that fixes the two issues you mentioned (an infinite function recursion issue and an issue with the "plussers" function for IQTensors). I came up with an even simpler implementation of plussers that no longer uses std::map, so that's nice. I added a unit test based on the code example you provided, so that ought to prevent the bug from recurring in the future.

Of course let me know if there's still an issue or if you spot another problem.

Miles