Possible bug in on-site terms - AutoMPO (C++)

Question

I am noticing a bug when I use AutoMPO for building the following bosonic Hamiltonian
$$
H = \sum_j n_j n_j
$$
Essentially, it looks like it does not do the square correctly. For example, if I run the following code (I don't impose conservation of particles, since in the actual model I have been investigating it is not conserved)

int size = 2;
int cutoff = 3;
SiteSet sites = Boson(size,{"ConserveQNs",false,"MaxOcc=",cutoff}); 
auto ampo = AutoMPO(sites); 
for(int j = 1 ; j <= size ; j++) ampo += 1 , "N", j , "N", j ;
MPO H = toMPO(ampo,{"Exact=",true});
PrintData(H(1));

The output does not appear to be n^2:

H(1) =
ITensor ord=3: (dim=4|id=952|"n=1,Site,Boson") (dim=4|id=952|"n=1,Site,Boson")' (dim=3|id=827|"l=1,Link")
{norm=8.60 (Dense Real)}
(1,1,1) 1.0000000
(2,2,1) 1.0000000
(3,3,1) 1.0000000
(4,4,1) 1.0000000
(2,2,2) 2.0000000
(3,3,2) 4.0000000
(4,4,2) 6.0000000
(2,2,3) 1.0000000
(3,3,3) 2.0000000
(4,4,3) 3.0000000

Am I missing something?

Thank you in advance!

Comments

Thanks and good catch: I think it is a bug. We have handled this case properly in the Julia version of ITensor but the C++ AutoMPO code is quite a bit more complicated and slightly older. I was able to reproduce this and will look into the best way to fix it.

Answer 1

Hi, so I checked some more on this. I think the code is working properly. It's hard (actually, impossible) to see from just looking at the first MPO tensor whether the MPO is correct because what operator is made by the MPO depends on all the MPO tensors.

Here is a check you can do to tell whether in this case the MPO returned equals the requested operator:

auto N2_1 = mapPrime( op(sites,"N",1)*prime(op(sites,"N",1)) ,2,1);
auto N2_2 = mapPrime( op(sites,"N",2)*prime(op(sites,"N",2)) ,2,1);
auto Id_1 = op(sites,"Id",1);
auto Id_2 = op(sites,"Id",2);

printfln("Norm difference = %.12f",norm(H(1)*H(2)-N2_1*Id_2-Id_1*N2_2));

The first two lines make the operator N^2 on each site. Then the last line subtracts (N^21 + N^22) from the operator defined by the MPO, which can be obtained by contracting all the MPO tensors together. When I run a code with your AutoMPO input at the top and this check at the end I find that the Norm difference as printed above is zero.

Please let me know if that doesn't answer your question or if you have a follow-up question.

Thanks,
Miles