dmrg function not accepting MPO obtained by using nmultMPO

Question

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?

Answer 1

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

Comments

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

---

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

---

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.

---

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.

---

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

---

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

---

Hi Prakash, sorry about the slow reply.

So DMRG in ITensor accepts an MPO where each MPO tensor has one primed (prime-level 1) and one unprimed (prime-level 0) physical index. In the ITensor library, as you may know, there is no ordering of the tensor indices which is ever required or fixed. So in that sense, for the purposes of DMRG, there's no such thing as an HH_(a'a) MPO versus an HH_(aa') MPO: they are the same thing.

(Of course if you print ITensors you will see the indices in a certain order. This is because we have to have some internal rule of how the indices are ordered and this is reflected in the printing. But in principle, ITensors could change their index ordering at any time, and as long as the elements inside are also permuted to match this new ordering, all of the other ITensor code should still work and give exactly the same result.)

Thus, if you are getting an error or the wrong answer, then it's not coming from the ordering of the indices of your MPO, but rather either because the virtual/link indices don't match each other, the physical indices don't match each other or those of the MPS, or because the numbers inside the MPO are wrong. (Or the DMRG calculation could just not be converging.)

Best,
Miles