Spinless Fermionic Correlator Across a 1-d Lattice

Question

I'm failing to generalize the example in the sample code Measure a correlator from an MPS wavefunction to assign the entries across the lattice in the C++ version. When applying to say, the Kitaev wire model, I can create the MPS and calculate one correlation, but when I try to generalize to measuring it across the lattice and assigning that to a 1-d array, I get an error Wrong number of IndexVals passed to elt/eltC. (expected 0, got 4), with that error having been expected 2 at some version for debugging. Additionally, when I attempt to do this, I get the following before that error message:
inds () = (dim=2|id=943|"n=2,Site,Fermion") <Out> 1: 1 QN({"Pf"}, 0, -2}) 2: 1 QN({"Pf"}, 1, -2}) (dim=2|id=173|"n=2,Site,Fermion") <In> 1: 1 QN({"Pf"}, 0, -2}) 2: 1 QN({"Pf"}, 1, -2}) (dim=2|id=173|"n=2,Site,Fermion") <Out> 1: 1 QN({"Pf"}, 0, -2}) 2: 1 QN({"Pf"}, 1, -2}) (dim=2|id=943|"n=2,Site,Fermion") <In> 1: 1 QN({"Pf"}, 0, -2}) 2: 1 QN({"Pf"}, 1, -2})
Which clearly has something to do with the QN Objects, but I'm not sure exactly what. It doesn't help that I haven't used C/C++ in a long time, so having most everything listed as auto for the datatype is somewhat dissatisfying. I attempted to use typeid to determine what I should be using for datatypes in these examples, but that failed fantastically. My guess is that there is an issue when I attempt to put the left and right index linking in multiple times? Then the operators that I stick into elt are defined properly?

I've attached two basic codes, a "Good" single correlation calculation, and a "Bad" multiple-site calculation. To keep things simple, I haven't tried to optimize the bad example, which will come after I get something to work. You'll also notice that I assign the correlation for site jj=3 to element 3 of the 1D array but L=20 here so it's not something I'm worried about. Though it is so embarrassing when my Fortran slips out in public.

Here is the "Good" implementation that runs, computes the ground state energy, and spits out a single correlation:

`
#include "itensor/all.h"
using namespace itensor;

   int main()
   {

    int L = 20;
       auto psites = Fermion(L,{"ConserveNf",false});
    auto mu = 0.1;
    auto t = 1.0 ;
       auto Delta = 1.0;

          // Create the operators for building the Hamiltonian
       auto ampo = AutoMPO(psites);
    for( auto ii:range1(L-1))
       {
           ampo += mu, "N", ii;
                 // ampo += mu, "N"-0.5, ii;    
        ampo += -t, "Cdag", ii, "C", ii+1; 
        ampo += -t, "Cdag", ii+1, "C", ii;

        ampo += Delta, "Cdag", ii, "Cdag", ii+1; 
        ampo += Delta, "C", ii+1, "C", ii;
    }

          // Add the final site's chemical potential term
    ampo += mu,"N",L;


 // Generate the Hamiltonian
    auto H = toMPO(ampo);
    auto state = InitState(psites);

   // Set the initial state to be a product state
   // such that it is within the symmetry sector 
   // OLD for Fermion Hopping: auto state = InitState(psites,"Emp");
for(int is = 1; is <= L; ++is)
{
    state.set(is, is%2==1 ? "1" : "0");
}
auto psi0 = randomMPS(state);

for( auto ii : range1(L) ) if( ii%2==0 ) state.set(ii,"Occ");
   // Set the number of DMRG sweeps
auto sweeps = Sweeps(5); 
   // Set the maximum bond dimension per sweep
   // Setting it low initially allows the computation to be much faster
sweeps.maxdim() = 10,20,100,100,200; 
   // Set the cutoff for the bond dimension
sweeps.cutoff() = 1E-10;

   // Run the DMRG
auto [energy,psi] = dmrg(H,psi0,sweeps,"Silent");
   // Calculate the energy and state

println("Ground state energy is ", energy) ;

auto ii =  1;
auto jj = 3;
auto Adag_ii = op(psites,"Adag",ii);
auto A_jj = op(psites,"A",jj);

   // 'gauge' the MPS to site i
   // any 'position' between i and j, inclusive, would work here
psi.position(ii);

auto psidag = dag(psi);
psidag.prime();

     // index linking ii to ii-1:
     auto lii_1 = leftLinkIndex(psi,ii);
auto Cij = prime(psi(ii),lii_1)*Adag_ii*psidag(ii);
for(int kk = ii+1; kk < jj; ++kk)
{   
    Cij *= psi(kk);
    Cij *= op(psites,"F",kk); // J-W string - need to add the "F" argument in between
    Cij *= psidag(kk);
}
// index linking jj to jj+1
auto ljj = rightLinkIndex(psi,jj);
Cij *= prime(psi(jj),ljj);
Cij *= A_jj;
Cij *= psidag(jj);

float GFjjresult    = elt(Cij); // use eltC if expecting complex

std::cout << GFjjresult << ' ';

   return 0; 
   }

`

And here is the "Bad" implementation where everything goes to heck:
`
#include "itensor/all.h"
using namespace itensor;

   int main()
   {

   int L = 20;
   auto psites = Fermion(L,{"ConserveNf",false});
   auto mu = 0.1;
   auto t = 1.0 ;
   auto Delta = 1.0;

   // Create the operators for building the Hamiltonian
   auto ampo = AutoMPO(psites);
   for( auto ii:range1(L-1))
{
       ampo += mu, "N", ii;
          // ampo += mu, "N"-0.5, ii;   
    ampo += -t, "Cdag", ii, "C", ii+1; 
    ampo += -t, "Cdag", ii+1, "C", ii;
   //   
    ampo += Delta, "Cdag", ii, "Cdag", ii+1; 
    ampo += Delta, "C", ii+1, "C", ii;
}

   // Add the final site's chemical potential term
   ampo += mu,"N",L;


   // Generate the Hamiltonian
   auto H = toMPO(ampo);
   auto state = InitState(psites);

   // Set the initial state to be a product state
   // such that it is within the symmetry sector 
   // OLD for Fermion Hopping: auto state = InitState(psites,"Emp");
   for(int is = 1; is <= L; ++is)
   {       
           state.set(is, is%2==1 ? "1" : "0");
   }
   auto psi0 = randomMPS(state);

   for( auto ii : range1(L) ) if( ii%2==0 ) state.set(ii,"Occ");
   // Set the number of DMRG sweeps
   auto sweeps = Sweeps(5); 
   // Set the maximum bond dimension per sweep
   // Setting it low initially allows the computation to be much faster
   sweeps.maxdim() = 10,20,100,100,200; 
   // Set the cutoff for the bond dimension
   sweeps.cutoff() = 1E-10;

   // Run the DMRG
   auto [energy,psi] = dmrg(H,psi0,sweeps,"Silent");
   // Calculate the energy and state

   println("Ground state energy is ", energy) ;

   auto ii =  1;
   Real GF1jj[L];
   auto jj = 3;
auto Adag_ii = op(psites,"Adag",ii);
auto A_jj = op(psites,"A",jj); 

   // 'gauge' the MPS to site i
   // any 'position' between i and j, inclusive, would work here
psi.position(ii);

auto psidag = dag(psi);
psidag.prime();

     // index linking ii to ii-1:

     for(int jj=ii+1; jj < L; ++jj)
    {
            auto lii_1 = leftLinkIndex(psi,ii);
    auto Cij = prime(psi(ii),lii_1)*Adag_ii*psidag(ii);
    for(int kk = ii+1; kk < jj; ++kk)
    {   
        Cij *= psi(kk);
        Cij *= op(psites,"F",kk); // J-W string - need to add the "F" argument in between
        Cij *= psidag(kk);
       }
    // index linking jj to jj+1
    auto ljj = rightLinkIndex(psi,jj);
    Cij *= prime(psi(jj),ljj);
    Cij *= A_jj;
Cij *= psidag(jj);

Real GFjjresult = elt(Cij); // use eltC if expecting complex


   GF1jj[jj] = GFjjresult;
}

    std::cout << GF1jj[3] << ' ';


   return 0; 

   }

`

Comments

I attempted to add ``Cij.noPrime();`` or ``Cij=noPrime(Cij);`` at the end of the loop to resent any priming/summation issues, but that did not work (tried both, but one at a time). Also tried writing out a ``jj=ii`` case and not having variable types declared inside the for loop. Neither seemed to fix this issue.

---

Also attempted to noPrime  both psi and psidag at the end of the loop iteration to no avail. It's clearly not constructing Cij correctly since that's what I feed into elt, but I don't see understand how it's constructing Cij with too many IndexVals if I unprime and reprime with the new desired sites.

---

The issue seems to be the definition of the operator A_jj with the dummy index jj defined before the loop, but not redefined as jj is changed, so when Cij was stuck into elt, it was poorly defined. Will post a working and better commented code once I have it working. No promises on it being optimized, just working.

---

Hi Jared, thanks for the question. I just got a chance to take a look at it, and yes that's the issue I found (that A_jj was defined to act on site 3, whereas in the code at some point you grab the MPS tensors for site 2 and try to compute an expectation value with them around A_jj, but for a different jj, 2 versus 3).

In general, my comments to you would be to:

(1) I don't recommend reusing the same variable for multiple purposes. I would personally use a loop counter with a name like "n" and then define jj differently, like n plus some constant or else give jj a starting value and then increment jj at the end of each turn of the loop, but not tie it to the loop counter.

(2) I would recommend printing out a lot more things (maybe you were doing this in your own personal version of the code). What I mean is putting in lines like:
Print(psi(jj));
Print(A_jj);

etc. Note to obtain this "Print" macro you will also need the line
#include "itensor/util/print_macro.h"
at the very top section of your code.

Printing everything out would have immediately showed you that A_jj had the 3rd site index whereas psi(jj) has the second site index.

Hope that helps!

---

Hi Miles,

Thank you for the comments. I certainly intended jj to be used for the same purpose (as a dummy for the site index for the annihilation operator), but clearly got a little sloppy when trying to extend the example. I have a working version now that I will post for anyone else who runs into a similar issue once I finish cleaning it up some.

I used print statements to isolate the issues, but did not know about the Print macro. It's always a bit frustrating to simultaneously learn the language and a package written by others.

Thanks again,
Jared

---

Hi Jared, yes I certainly understand that if you’re learning C++ then the ITensor specific features are hard to find out about and to learn at the same time.

Just to follow up - you got the code working the way you wanted now? If so I’ll write an answer below to mark the question as answered. If you did post the code that would certainly be welcome.

Thanks,
Miles

---

Hi, Miles,

Yes, it seems to be working. I will a couple of tests before I post it. My goal is to post code that explains the idea of what's going on, so I want to go through it a couple more times before posting it.

Thanks,

Jared

---

That sounds good, thank you.

Answer 1

(Please see the discussion above.)