How to get the MPO of H^2 from the one of H ?

Question

Hi everyone,

I am trying to get the MPO representation of H^2 from the one of H using the
nmultMPO function but I am unable to get the correct result. (Notice H is long ranged Hamiltonian) Below is the piece of code I am using:

auto sites = SpinHalf(N,{"ConserveQNs=",false}); //make a chain of N spin 1/2's
auto ampo = AutoMPO(sites);

 for(auto i : range1(N-1))
 {
        for(auto j=i+1;j<=N;j++)
        {
            ampo += (i+j),"S+",i,"S-",j;
            ampo += (i+j),"S-",i,"S+",j;
        }
 } 
auto H = toMPO(ampo); 
auto res= nmultMPO(H, H,{"Cutoff",1e-8}); 

I get the following error:
libc++abi.dylib: terminating with uncaught exception of type itensor::ITError: doTask not defined for task Mult and storage type Combiner

Your help will be greatly appreciated.

Regards,
Yas

Comments

How about writing down the form of H^2 and constructing it directly using AutoMPO?

---

Thank you for your comment hermit0308.
Well, I would prefer not as I want to then compute successive power of H.

Answer 1

Hi Yas,

The new version of nmultMPO in V3 works a little bit differently from the one from V2. Now, the input is a little bit more "literal", in that it finds the indices common to the two input MPOs, contracts over those, and returns an MPO with the indices that are not common between the two MPOs.

In short, you should call the function like this:

auto res = nmultMPO(H, prime(H),{"Cutoff",1e-8});
res.mapPrime(2,1);

to get the same thing as ITensor V2.

The problem was that the function wasn't working since it was automatically contracting over all of the site indices of the input MPO. However, this should have a better error message, I will add a check for that.

To explain why the behavior changed, the reason is that there are many MPOs you could make (such as in more general tensor networks) that don't have pairs of indices (s1,s1',s2,s2',...), and that wasn't allowed in the previous version. The new behavior makes the code a bit more general (and I personally feel is a bit more intuitive, since now contracting MPOs works more like contracting ITensors).

Cheers,
Matt

Comments

Hi Matt,

Many thanks for your answer: it is now working.

Cheers,
Yas