# dmrg function not accepting MPO obtained by using nmultMPO

+1 vote

Hi,
I am using itensor v3 to compute ground state of (H-E)^2=H^2-2EH+E^2, where H is Fermi Hubbard Hamiltonian as defined in http://itensor.org/docs.cgi?page=classes/autompo&vers=cppv3 except I am calculating it in 2D square lattice and E is some constant. In order to construct H^2 I used nmultMPO function flowing way:
auto H = toMPO(ampo);
auto Hn = toMPO(ampo);
auto HH = MPO{};
nmultMPO(prime(H),Hn,HH,{"MaxDim",1000,"Cutoff",1E-14});
so constructed HH (i.e. H^2) works fine with operation like inner(psi,HH,psi) but calling it upon dmrg function returns following message:
dim(inds(phi.front())) = 16
A.size() = 64
From line 137, file /data/pn1_1/psharma/ITensor-v3/itensor/iterativesolvers.h

davidson: size of initial vector should match linear matrix size

davidson: size of initial vector should match linear matrix size
Aborted (core dumped)
while dmrg function works perfectly fine with H and same initial MPS (used for HH).
Is there a way to fix this issue or do there exist any other simple way to find ground state of my problem?

Thanks for the question. So the issue here is that the site indices of the MPO HH do not have the prime levels expected for an MPO that is used as input to the dmrg function. If you print out HH, you'll see that each tensor has a site index with prime level 0 and a site index with prime level 2. That's in contrast to a usual MPO for a Hamiltonian in ITensor, where the prime levels of the site indices are 0 and 1.

The reason it's 0 and 2 in your case is that prime(H) raises the prime levels from 0 and 1 to 1 and 2, then contracting that MPO with H contracts over the prime-level-1 indices, leaving only the prime level 0 and 2 indices.

So a solution is just to add the line

HH.mapPrime(2,1);


after you make HH, which will map any prime-level-2 indices to have prime level of 1 instead.

I'd encourage you to print out the MPO you are inputting into dmrg to make sure it has the structure and prime levels you want.

The more surprising thing here is that inner(psi,HH,psi) works. It might be a good idea in the future for us to make that more restrictive, but the current way that function works is that as long as one set of site indices of the MPO matches the MPS on the right, then the other site indices get automatically modified to match those of the MPS on the left. This is mentioned in the documentation here: http://itensor.org/docs.cgi?vers=cppv3&page=classes/mps_mpo_algs
But I agree it can be a bit surprising in cases like this.

Best,
Miles

commented by (160 points)
Hello Miles,
Thank you for your reply. Sorry that I forgot to mention in my previous message that I actually had implemented mapPrime function to lower the prime level. As you mentioned on your reply, it surprisingly works for inner(psi,HH,psi) even without lowering the  prime level to 1. Upon calling this HH MPO for dmrg calculation, no matter whether you lower the prime level or not, returns the same error message except the value of  A.size() (in the second line of error message). It is complaining about the size of the initial vector (which I guess is initial MPS) and that of linear matrix (which is HH MPO I think), which is I don't understand clearly why these sizes are compatible in the case of H but not in the case of HH obtained my multiplying two MPOs (H and Hn, note H=Hn).
Sorry for bothering you again but, can you please give me some idea what should be a better way to solve this issue?

best regards,
Prakash
commented by (70.1k points)
Hi Prakash,
Thanks for following up. So are you saying that adding a line HH.mapPrime(2,1); before passing HH to DMRG does not solve the issue? I made a sample code for myself and adding this line did solve the issue for me. So in order to know why it is not working for you, I would need to know all of the steps you are doing, and how they differ from mine. Could you therefore post the relevant part of your code (from making the MPO to calling DMRG, at least)?

Also I would recommend that you print out the MPO you are passing to the DMRG function to see if it has the expected prime level structure. Each tensor should have two site indices, with one site index having prime level 0 and the other having prime level of 1. Printing out the MPO should let you see whether this is the case.

Thanks,
Miles
commented by (160 points)
Hello Miles,
I have attached my entire code below. You said your code worked  after lowering the prime level; however, it's not working for me and I don't understand why. Its possible I might be doing something wrong but I don't see it. Please suggest me the correct solution.

#include "itensor/all.h"

using namespace itensor;

int main(int argc, char* argv[])
{

auto Nx = 3,
Ny = 3;
auto N = Nx*Ny;

auto sites = Electron(N);

auto t = atof(argv[1]);
auto U = atof(argv[2]);

auto ampo = AutoMPO(sites);
auto lattice = squareLattice(Nx,Ny,{"YPeriodic=",true});

printfln("square lattice NN pairs");

for(auto& bnd : lattice)
{
printfln("Bond from site %d -> %d",bnd.s1,bnd.s2);
}

for(auto j : lattice)
{
ampo += -t,"Cdagup",j.s1,"Cup",j.s2;
ampo += -t,"Cdagup",j.s2,"Cup",j.s1;
ampo += -t,"Cdagdn",j.s1,"Cdn",j.s2;
ampo += -t,"Cdagdn",j.s2,"Cdn",j.s1;
}
for(auto j : range1(N))
{
ampo += U,"Nupdn",j;
}
auto H = toMPO(ampo);
auto Hn = toMPO(ampo);

auto HH = MPO{};
nmultMPO(prime(H),Hn,HH,{"MaxDim",1000,"Cutoff",1E-14});
HH.mapPrime(2,1);

//Print(HH.A(1).inds());
//auto s3 = sites(3);
//Print(hasInds(HH(3),{s3,prime(s3,1)}));

auto state = InitState(sites);
for(auto j : range1(N)){state.set(j,(j%2==1 ? "Up" : "Dn"));  }
auto psi0 = MPS(state);

PrintData(inner(psi0,HH,psi0));  //for comparison
PrintData(inner(H,psi0,H,psi0));

auto sweeps = Sweeps(5);
sweeps.maxdim() = 10,20,100,100,200;
sweeps.noise() = 1E-7,1E-8,1E-10,0;
sweeps.cutoff() = 1E-6;

auto [energy,psi] = dmrg(HH,psi0,sweeps,{"Quiet=",true});

printfln("GS energy = ",energy);

return 0;
}

Thanks a lot.

-Prakash.
commented by (160 points)
Just now I found it works if I switch the position of prime(H) and Hn inside nmultMPO function i.e. nmultMPO(H,prime(Hn),HH) works but the one in my above code doesn't.
commented by (70.1k points)
Hi Prakash,
That makes sense. Did you draw a diagram of Hn and prime(H) and how nmultMPO is multiplying them together? Are you now getting the results you expect? I know these things can be tricky.

Miles
commented by (160 points)
Hello Miles,
I think we get HH_(a'a) in the first case and HH_(aa') in the later (here I am assuming a and a' indices after lowering the prime level). So, dmrg function accepts only MPO of type H_(aa')?
If I am thinking wrong please correct me.

Best,
Prakash
commented by (70.1k points)