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

Do you know why I had to use "removeQNs"? My sites conserve Sz.

---

Did you have to use removeQNs? I did not see this in the code above or in our discussion so far. Where did you have to use it? I don’t think I can give a more helpful answer until you let me know more context and more details of your question - thanks.

---

Here are the parts of the code that are relevant: https://pastebin.com/TGeDqwTh
Here is the full code: https://pastebin.com/6iVvr4Ht

As a summary, I create wavefunctions that are products of dimers (where the dimer is defined as all 4 states of the Sz=0 sector of 2 sites), and I store them in an MPS vector. I compute their energies for the Hubbard Hamiltonian which is stored in a MPO vector (vector because I use many U's).  This works fine. But now, I want to compute these matrix elements that I showed in the  question. In order to do that, if I keep conserving Sz, I have to removeQNs in the computeG function. Do you need any more info?

---

Something I forgot to mention is that I optimize the coefficients inside each dimer, which means that each cluster has the exact ground state (with a mean-field of the rest of the clusters). This means that I know that some matrix elements from the equation at the question are 0 and some others non-zero. The problem is that I am getting all of them zero...

---

Sorry, this is rather complicated code to go all through so I won't be able to without further discussion to focus on your specific question. When you say you had to use removeQNs, do you mean you did not have it before and you got a bug, then you added it and it worked? Or did you have a specific reason for adding it, such as mixing QN-conserving ITensors with dense ITensors?

---

So I have a QN conserving MPS and an MPO (the Hamiltonian) that is also QN conserving. The cc operator, for example,  is Cdagup,1,Cup,2. When I apply it to the QN-conserving MPS, I get no error, but when I compute inner(MPS, H, cc_MPS), I get a "Mismatched IQIndex arrows" error. This is without removeQNs. If I add removeQNs, it computes a number, but I do not think it is the correct number, so there is still a problem. Could there be a problem in the way I define ampo? Does  the order of the site indices matter? Should I use Electron instead of Hubbard?

---

I see - thanks for explaining that. So most likely what's happening is that when you modify one of the MPS, the arrows of the site indices are getting changed to the wrong direction. I would guess this is due to an improper construction of the operators you're using or perhaps how you are acting them on the MPS. We can discuss which way the arrows should go (all site indices Out for a ket MPS) and for operators (one In and unprimed, one Out with prime-level 1) in more detail and the reasons, but perhaps you already understand that part well?

So what I would mainly suggest is to print out the steps of when you apply the operators to the MPS and draw diagrams of what's happening at each step (I often do this with my own code). A good pattern to use is to put the macro PAUSE; after you print, which will pause the code there until you hit Enter so you can study the printout without having to scroll up through a lot of text, or so you can pause each step of a loop and inspect what is happening slowly.

I would treat this arrow issue as a very serious issue. Using removeQNs to avoid it is not really fixing the problem but just neglecting it and likely giving a wrong answer (because as you may know we don't enforce Index arrow rules for dense tensors, though in the future we will offer this optionally).

---

So here is the output that I get:
MPS:
ITensor ord=2:
(4|id=238|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=559|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})
{norm=2.00 (QDense Real)}

ITensor ord=2:
(4|id=559|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=233|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})
{norm=1.00 (QDense Real)}

ITensor ord=2:
(4|id=573|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=320|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})
{norm=2.00 (QDense Real)}

ITensor ord=2:
(4|id=320|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=258|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})
{norm=1.00 (QDense Real)}
-----------------
I think the MPS looks fine to me.
-------------------------------
MPO, which is the Cdagup, 1, Cup, 2:
ITensor ord=3:
(3|id=279|l=1,Link) <Out>
  1: 2 QN()
  2: 1 QN({"Nf",-1,-1},{"Sz",1})
(4|id=238|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=238|n=1,Elec,Site) <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})
{norm=2.45 (QDense Real)}

ITensor ord=4:
(3|id=279|l=1,Link) <In>
  1: 2 QN()
  2: 1 QN({"Nf",-1,-1},{"Sz",1})
(2|id=165|l=2,Link) <Out>
  1: 2 QN()
(4|id=233|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})
(4|id=233|n=2,Elec,Site) <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})
{norm=3.16 (QDense Real)}

ITensor ord=4:
(2|id=165|l=2,Link) <In>
  1: 2 QN()
(2|id=616|l=3,Link) <Out>
  1: 2 QN()
(4|id=573|n=3,Elec,Site) <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=573|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})
{norm=2.83 (QDense Real)}

ITensor ord=3:
(2|id=616|l=3,Link) <In>
  1: 2 QN()
(4|id=258|n=4,Elec,Site) <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=258|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})
{norm=2.00 (QDense Real)}
-----------------
What do you think of that?
--------------------------
And lastly, the MPO acting on the MPS:
ITensor ord=2:
(2|id=275|U,Link) <Out>
  1: 1 QN({"Nf",2,-1},{"Sz",0})
  2: 1 QN({"Nf",3,-1},{"Sz",1})
(4|id=238|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})
{norm=0.71 (QDense Real)}

ITensor ord=3:
(1|id=373) <Out>
  1: 1 QN({"Nf",2,-1},{"Sz",0})
(4|id=233|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=275|U,Link) <In>
  1: 1 QN({"Nf",2,-1},{"Sz",0})
  2: 1 QN({"Nf",3,-1},{"Sz",1})
{norm=1.41 (QDense Real)}

ITensor ord=3:
(4|id=985|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=573|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})
(1|id=373) <In>
  1: 1 QN({"Nf",2,-1},{"Sz",0})
{norm=1.00 (QDense Real)}

ITensor ord=2:
(4|id=258|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=985|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})
{norm=2.00 (QDense Real)}
----------------------
There is something in the 3rd iTensor that doesn't look right to me...


What do you think?

---

At a glance it looks ok. You'll need to remove the primes from the result MPS to restore it back to ITensor's standard MPS format (no primes on any indices) before doing a subsequent operation. But I don't see any arrow issues right away.

We do know that you got an arrow error, though, so you'll have to keep looking to find out which line of code gave that error. If you run your code in a debugger, for example, you can track back to the exact step.

---

For the first part, you just mean putting A.ref(j) = noPrime(A.ref(j));, right?

Now regarding the rest. There is another problem now... I tried "make debug" and the file-g crashes at a completely different step... It crashes when I evaluate <Psi|H|Psi> and the error is
Non-QN index = (1|id=0)
From line 547, file indexset.cc

IndexSet: cannot mix QN and non-QN Indices

How can that be the case if I only have sites that conserve the QNs...?

---

Yes, sort of, to the first part. Though it would be much better to call `A.noPrime();` to remove all of the primes from the MPS `A`. Apart from being shorter and not requring a loop, it would also preserve any orthogonality properties of `A`. (MPS objects in ITensor "remember" which tensors are left- and right-orthogonal to skip past certain operations and make your code run faster.)

I'm not sure about the other bug without more context. You'll need to print out the tensors prior to that line to see if they do carry QNs and whether everything is defined ok. It could even be that one of the tensors isn't defined at all (since I see an id=0 which using means a default-constructed Index).

---

After lots of trial error, I tried something and it works, but it gives the same result as removeQNs...

Instead of inner(Phi_0, cc, h_mps); I wrote inner(cc_mps,h_mps);, where cc_mps = applyMPO(cc,MPS); and h_mps=applyMPO(h,MPS);

Isn't this the same, though? Or the bras and the kets are handled differently?

---

It should give the same result, or close, though it uses a different algorithm which is both more computationally costly and less accurate in general.

I would definitely recommend not just doing trial and error, but instead printing out your tensors at each step and/or using a debugger to find the real reason for the arrow error.

---

I am not sure if I completely understand what you mean.

I have printed the MPS before the applyMPO() and after it. The error is inside inner(Phi_0,cc,h_mps);, so are you saying that I should find in which line of inner() the code  is crashing?

Another thing I can print is with PrintData()

Phi_0.ref(pnt) =
ITensor ord=2:
(4|id=248|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=644|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})
QDense Real {4 blocks; data size 4}
Block: 1 1
(1,1) 1.0000000
Block: 2 2
(2,2) 1.0000000
Block: 3 3
(3,3) 1.0000000
Block: 4 4
(4,4) 1.0000000

cc_mps.ref(pnt) =
ITensor ord=2:
(2|id=669|U,Link) <Out>
  1: 1 QN({"Nf",2,-1},{"Sz",0})
  2: 1 QN({"Nf",3,-1},{"Sz",1})
(4|id=248|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})
QDense Real {2 blocks; data size 2}
Block: 2 3
(2,3) 0.4351507
Block: 1 4
(1,4) -0.557343

Phi_0.ref(pnt) =
ITensor ord=2:
(4|id=644|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=135|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})
QDense Real {4 blocks; data size 4}
Block: 4 1
(4,1) 0.4351563
Block: 3 2
(3,2) -0.557361
Block: 2 3
(2,3) 0.5573434
Block: 1 4
(1,4) 0.4351507

cc_mps.ref(pnt) =
ITensor ord=3:
(1|id=207) <Out>
  1: 1 QN({"Nf",2,-1},{"Sz",0})
(4|id=135|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=669|U,Link) <In>
  1: 1 QN({"Nf",2,-1},{"Sz",0})
  2: 1 QN({"Nf",3,-1},{"Sz",1})
QDense Real {2 blocks; data size 2}
Block: 1 1 1
(1,1,1) 1.0000000
Block: 1 2 2
(1,2,2) 1.0000000

Phi_0.ref(pnt) =
ITensor ord=2:
(4|id=820|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=476|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})
QDense Real {4 blocks; data size 4}
Block: 1 1
(1,1) 1.0000000
Block: 2 2
(2,2) 1.0000000
Block: 3 3
(3,3) 1.0000000
Block: 4 4
(4,4) 1.0000000

cc_mps.ref(pnt) =
ITensor ord=3:
(4|id=590|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=820|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})
(1|id=207) <In>
  1: 1 QN({"Nf",2,-1},{"Sz",0})
QDense Real {4 blocks; data size 4}
Block: 4 1 1
(4,1,1) -0.435141
Block: 3 2 1
(3,2,1) -0.557320
Block: 2 3 1
(2,3,1) 0.5573798
Block: 1 4 1
(1,4,1) -0.435171

Phi_0.ref(pnt) =
ITensor ord=2:
(4|id=476|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=616|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})
QDense Real {4 blocks; data size 4}
Block: 4 1
(4,1) 0.4351714
Block: 3 2
(3,2) -0.557380
Block: 2 3
(2,3) 0.5573202
Block: 1 4
(1,4) 0.4351411

cc_mps.ref(pnt) =
ITensor ord=2:
(4|id=616|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=590|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})
QDense Real {4 blocks; data size 4}
Block: 1 1
(1,1) 1.0000000
Block: 2 2
(2,2) 1.0000000
Block: 3 3
(3,3) 1.0000000
Block: 4 4
(4,4) 1.0000000

Phi_0 is the MPS that I have and cc_mps is cc acting on Phi_0. Should I look for something specifically?

---

Hi, yes so no need to do PrintData unless you think the problem has something to do with the entries of one or more tensors (e.g. one of the tensors is undefined or zero, say). But that's not likely the case here. So the printing (just using Print) is to look at the indices, to see if they match (same id's, prime levels, and tags) and have compatible arrows (in arrows are only allowed to contract with out arrows).

For the call to inner(Phi_0,cc,h_mps), if that's where the crash is coming from then there are only two main possibilities:
1. there is a bug in the inner function
2. the inputs are not of the form that `inner` expects, for some reason

#2 is the most likely, so the best thing to do next is to print out Phi_0, cc, and h_mps and look at the indices of their tensors. Also after printing them out, draw the diagram of the expectation value <Phi_0|cc|h_mps> which is what inner computes to see if the indices can all be connected in the way expected by the inner function. Note that inner takes the Hermitian conjugate of the first argument (Phi_0 here) for you, so Phi_0 and h_mps should both be well-formed MPS with the same site indices as each other and as h_mps, and these site indices should all have Out arrows with prime-level zero (no prime). Is that the case?

If that's all the case and you are still finding a crash in inner, then I'll need a minimal working example from you to see if I can reproduce the bug in inner so I can fix it.

Hope that helps!

Miles

---

Thanks again so much for your help!

Here I have a minimal code: https://pastebin.com/Z4EucnLz
I have also put comments. "State" contains the state Phi_0 and "H" contains the hubbard Hamiltonian. The rest is in the code.

You do not really need to look the rest of the functions, only computeG which is the function that computes the matrix elements of the question.

---

I tried the following:

                                        MPS h_mps = applyMPO(h,Phi_0); MPS cc_mps = applyMPO(cc,Phi_0);
                                        cc_mps.noPrime(); h_mps.noPrime();

                                        MPS h_cc_mps = applyMPO(h, cc_mps); MPS cc_h_mps = applyMPO(cc,h_mps);
                                        h_cc_mps.noPrime(); cc_h_mps.noPrime();

                                        double elem = -inner(Phi_0,h_cc_mps) + inner(Phi_0,cc_h_mps);

but I do not think that it gives the correct results, either. Because applyMPO turns both MPS's to prime indices instead of no-prime I guess I have to put the noPrime() after, right?

Any input on why the results seem wrong, though?

------------------------------------
But now I tried something different and the results make more sense. Could it be that I need to use Adag and A instead of Cdag and C?

---

Hi, I just took a look at the sample code. Unfortunately I'm not able to compile it yet, because I don't have the gsl library. Is that easy to install?

What is the status of the code at this point? Are you still getting errors / crashes? Or is it that you are not getting the result you expected?

When you say you use Adag and A instead of Cdag and C, do you mean inside AutoMPO or outside? The way AutoMPO handles fermions is different from the rest of the code (automatic inside AutoMPO; outside you must include Jordan-Wigner string). It can be a point of confusion but right now is the best we can do until we implement our new fermion system.

Miles

---

I dropped the dependency on GSL: https://pastebin.com/1SWpuCGY
Now you should be able to compile it easily.

If we do step by step applyMPO and after every applyMPO we do noPrime(), the code will work. But the results are not the ones I expect... More specifically, for the case of c!_1_up * c_2_up, this matrix element with the commutator should give 0, and for the case c!_2_up * c_3_up, this matrix element should not be zero (this is a special case for the wavefunction that I have).

I used Adag and A to define c!c, but I use Cdag and C to define the Hamiltonian. Anyway, I do no think it is correct, because I ran a few more tests and the answers still do not make lots of sense...

Once again, thank you very much!

---

Hi great. So this helped me to find one issue in the code, which is not quite a bug in ITensor, but more of a limitation. The inner function is requiring that the MPS (and MPO) you pass into it has to have a link index connecting each of the MPS tensors. I notice that State, being a product of dimers, does not have any link indices on the even-numbered bonds. So you will need to put these in and then that part of the code should work without any crashes.

(A convenient way to add dimension-1 Indices is just to make them like this:
auto i = Index(QN(),1,"i");
and then "stick" them onto a tensor T like this:
T *= setElt(i(1));
You'll have to stick dag(i) on the other MPS tensor.
)

In terms of why you are not getting the answer you expect, there could be many reasons, so it's hard to know why without more information. Could you possibly tell me a specific MPO and matrix element you are trying to compute, how each is defined precisely, then what answer you expect?

---

Do you mean adding something like that?

for (int ii=2; ii<=4; ii+=2) {
        auto ind = Index(QN(),ii,"i");
        Phi_0.ref(ii) *= setElt(ind(ii));
        Phi_0.ref(ii-1) *= setElt(dag(ind(ii)));
}

I tried it, but I am getting a "In findIndex: more than one Index found, consider using findInds instead" error...

---

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