Commutators of MPO's

Question

Hello!

I want to compute the following matrix elements:

The Hamiltonian (H) is already stored as an MPO and Phi_0 is also stored as an MPS. Below is the code I am trying right now:

        for (int i=0; i<N; i++)
    {
            for (int j=0; j<N; j++)
            {
               auto ampo = AutoMPO(sites);
               ampo += "Cdagup", i, "Cup", j;
               //ampo += "Cdagdn", i, "Cdn", j;
               MPO op = toMPO(ampo);

               auto C = nmultMPO(H[point],op,{"MaxDim",500,"Cutoff",1E-8}); cout << "os edo kala\n";
               auto D = nmultMPO(op, H[point], {"MaxDim",500,"Cutoff",1E-8}); cout << "os edo kala\n";
               gsl_matrix_set(G, i,j,-inner(states.at(ind),C,states.at(ind))+inner(states.at(ind),D,states.at(ind)));
            }
    }

I am getting the following error: Unassigned site in SiteStore

Any idea about why? Is it because I only have two indices in the C,D MPO's?

What I am trying to do is to basically implement orbital optimization for the Hubbard Hamiltonian. I do not know if there is another way to do that in iTensor.

Answer 1

Hi, so I'd recommend the following approach instead:

1. expand the commutator into two terms (H on the left, H on the right) and compute these separately and/or use a symmetry to relate one to the other (on paper)

2. for each of these terms, say <Phi|H c^\dagp cq |Phi> compute an MPS "ccPhi" by acting cq and c^\dagp onto Phi. You can do this by acting those operators on the site index of the appropriate MPS tensor and re-setting the prime level back to zero for those indices. You'll also need to act with Jordan-Wigner string on the sites between p and q. (Same kind of operation, but just act with the operator named "F".)

3. compute the overlap <Phi|H|ccPhi> by using the inner(Phi,H,ccPhi) function

That approach ought to be a lot faster than making an MPO for each c^\dagp cq pair. The only downside is needing to explicitly deal with Jordan-Wigner string though.

---

On the other hand, if you do prefer to stick with the MPO based approach, and if it's fast enough for you (at a minimum it's a good check on the other approach) then I'd recommend not using nmultMPO to multiply MPOs together, but instead take advantage of the fact that you are ultimately taking an expectation value with MPS. So you will find it works better to do applyMPO to multiply each MPO onto |Phi> until only one MPO is left, then use the inner function to compute the matrix element of that leftover MPO.

So act the "Cdagup",i,"Cup",j MPO onto |Phi> using applyMPO, to get |ccPhi>, then compute <Phi|H|ccPhi>.

Hope that works well for you!

Miles

Comments

Yes - just like that. However I think your for loop needs to start at 2 (and Id recommend making the ending condition based on the length of Psi_0 of course).  Did you print out the MPS afterward to make sure the tensors have the right structure?

---

The final structure is that:


ITensor ord=3:
(4|id=113|n=1,Elec,Site) <Out>
  1: 1 QN({"Nf",0,-1},{"Sz",0})
  2: 1 QN({"Nf",1,-1},{"Sz",1})
  3: 1 QN({"Nf",1,-1},{"Sz",-1})
  4: 1 QN({"Nf",2,-1},{"Sz",0})
(4|id=81|U,Link) <In>
  1: 1 QN({"Nf",0,-1},{"Sz",0})
  2: 1 QN({"Nf",1,-1},{"Sz",1})
  3: 1 QN({"Nf",1,-1},{"Sz",-1})
  4: 1 QN({"Nf",2,-1},{"Sz",0})
(2|id=197|i) <In>
  1: 2 QN()
{norm=2.00 (QDense Real)}

ITensor ord=3:
(4|id=81|U,Link) <Out>
  1: 1 QN({"Nf",0,-1},{"Sz",0})
  2: 1 QN({"Nf",1,-1},{"Sz",1})
  3: 1 QN({"Nf",1,-1},{"Sz",-1})
  4: 1 QN({"Nf",2,-1},{"Sz",0})
(4|id=276|n=2,Elec,Site) <Out>
  1: 1 QN({"Nf",0,-1},{"Sz",0})
  2: 1 QN({"Nf",1,-1},{"Sz",1})
  3: 1 QN({"Nf",1,-1},{"Sz",-1})
  4: 1 QN({"Nf",2,-1},{"Sz",0})
(2|id=197|i) <Out>
  1: 2 QN()
{norm=1.00 (QDense Real)}

ITensor ord=3:
(4|id=513|n=3,Elec,Site) <Out>
  1: 1 QN({"Nf",0,-1},{"Sz",0})
  2: 1 QN({"Nf",1,-1},{"Sz",1})
  3: 1 QN({"Nf",1,-1},{"Sz",-1})
  4: 1 QN({"Nf",2,-1},{"Sz",0})
(4|id=921|U,Link) <In>
  1: 1 QN({"Nf",0,-1},{"Sz",0})
  2: 1 QN({"Nf",1,-1},{"Sz",1})
  3: 1 QN({"Nf",1,-1},{"Sz",-1})
  4: 1 QN({"Nf",2,-1},{"Sz",0})
(4|id=627|i) <In>
  1: 4 QN()
{norm=2.00 (QDense Real)}

ITensor ord=3:
(4|id=921|U,Link) <Out>
  1: 1 QN({"Nf",0,-1},{"Sz",0})
  2: 1 QN({"Nf",1,-1},{"Sz",1})
  3: 1 QN({"Nf",1,-1},{"Sz",-1})
  4: 1 QN({"Nf",2,-1},{"Sz",0})
(4|id=856|n=4,Elec,Site) <Out>
  1: 1 QN({"Nf",0,-1},{"Sz",0})
  2: 1 QN({"Nf",1,-1},{"Sz",1})
  3: 1 QN({"Nf",1,-1},{"Sz",-1})
  4: 1 QN({"Nf",2,-1},{"Sz",0})
(4|id=627|i) <Out>
  1: 4 QN()
{norm=1.00 (QDense Real)}


It looks fine to me, isn't it? The even sites have "out" indices and the odd ones have "in" indices, and all of them are "i".

---

I see. So the first tensor has 3 indices which is the incorrect structure for a finite MPS. That is what I was thinking would happen since the loop starts at 1.

---

Yeah, I did it now, and it gives the same results as doing another applyMPO and the again noPrime. So let me explain now why I think I am getting wrong results.

What I am trying to achieve is to do orbital optimization. The formula that I have given is for computing the gradient for the orbital optimization.

In my case, the coefficients of the dimers are optimized (exact for each dimer) and so orbital rotations inside each dimer, meaning for ie. q=1 p=2, q=2 p=1,q=3 p=4, etc, this matrix element should be zero and for ie. q=1  p=3, etc we should get a value close to 1.  

This is not the case in my code... I am getting smalle values for nearly all p's and q's. I do not know the exact values, only roughly what I just mentioned.

The final code is this one: https://pastebin.com/tzZ4XG4R
In the code, for i=1, j=2, I should get zero (both spin up and spin down), and for i=1, j=3, I should get something large, as I previously said. The difference is that in the code I denote p and q as i and j respectively.

---

Hi yes, so it sounds like your code still has a flaw or bug in it. Do you think it is a bug in ITensor though? The last issue we discussed was definitely about a limitation of the ITensor `inner` function (requiring all of the links of an MPS to be there to work). But this issue sounds like it might be in the mathematics that your code is carrying out.

So in that case, all I can say to help is to do the following:
- print out key quantities to see if they have the values you expect
- test on small systems or in exact cases to see if the results match your expectations
- test small "components" of your code to see if they work correctly individually

Those are the steps I would take when writing a similar code. If you do suspect a part of ITensor is working improperly, though, or you aren't sure how to use it correctly please feel free to ask about that and/or post a new question.

Best,
Miles