Difficulty adding fermionic IQTensors with sum(psi1,psi2)

Question

I am attempting to create a trial wavefunction that is a superposition of two wavefunctions.

My code does roughly the following:

Hubbard sites;
auto sites1,sites2 = sites; // N is even
sites = Hubbard(N);
sites1 = sites;
sites2 = sites;

IQMPS psi1, psi2;

auto state1 = InitState(sites1);
auto state2 = InitState(sites2);
loop over sites;
      for odd sites, state1.set(i, "UpDn")
      for even sites, state2.set(i, "UpDn")
      other sites are empty

psi1 = IQMPS(state1);
psi2 = IQMPS(state2);
psi = sum(psi1,psi2);        // I also tried out sum(psi1, psi1) as a check

psi.norm();

The electron densities that result for psi are
8 8 8 8 8 8
whereas I was expecting
1 1 1 1 1 1

Notes:
The above code works fine for a single wavefunction. That is, if I just set psi = psi1 or psi = psi2, just the odd or just the even sites are filled, the electron density is
2 0 2 0 2 0
or
0 2 0 2 0 2

When I set psi = sum(psi1, psi1) the electron density is
8 0 8 0 8 0

The only portion of the code that isn't working, from what I can tell, is the sum( , ) function. Is there a bug in sum() that is not conserving quantum numbers? Is there another function that I can use to conserve quantum numbers while adding wavefunctions?

Comments

I forgot to mention my method for measuring the electron density. I used the method in the ITensor Samples folder, exthubbard.cc:

    loop over j
    psi.position(j);
    upd(j) = (dag(prime(psi.A(j),Site))*sites.op("Nup",j)*psi.A(j)).real();
    dnd(j) = (dag(prime(psi.A(j),Site))*sites.op("Ndn",j)*psi.A(j)).real();

and then output the sum of upd and dnd for each j.

After I do the git pull I will post whether this works or not.

Answer 1

Hi Jon,

I get the correct density (222222) with the latest version of ITensor. My code is

#include "itensor/all.h"

using namespace itensor;

int main()
{
auto sites = Hubbard(8);
auto sites1 = sites;
auto sites2 = sites;

IQMPS psi1, psi2;

auto state1 = InitState(sites1);
auto state2 = InitState(sites2);
for(int j = 1; j <= 8; ++j)
{
    if(j%2 == 1)
    {
        state1.set(j, "UpDn");
        state2.set(j, "Emp");
    }
    else
    {
        state2.set(j, "UpDn");
        state1.set(j, "Emp");
    }
}

psi1 = IQMPS(state1);
psi2 = IQMPS(state2);
auto psi = sum(psi1,psi2);        // I also tried out sum(psi1, psi1) as a check
Print(psi);
for(int j = 1; j <= 8; ++j)
{
    auto Ni = sites1.op("Ntot", j);
    auto ket = psi.A(j);
    auto bra = dag(prime(ket, Site));
    std::cout << (bra*Ni*ket).real() << std::endl;
}
return 0;
}

The last part gives electron density 2 on every site. Is that what you expect? Thanks.

Jin

Comments

Hi Jin,
Thanks a lot for posting a quality reply. I did notice that in your measurement loop, though, there is no call to psi.position(j) to ensure that the gauge/orthogonality center is at site j before doing the measurement. Can you check that if you include this line in your test code that it does continue to give the right answer?

---

Thanks Jin and Miles! I'll do a fresh git pull to use the latest revision.

---

Ok guys I got it working, thanks a lot!

---

Great! Ok thanks again Jin.

---

I fixed my code as follows.

I have an if() statement that enables "sites" to be read in from a file if such a file exists. If it does not exist, I create "sites" by setting it equal to "Hubbard(N_sites)". Defining "sites1" and "sites2" as type Hubbard before this if statement fixed the problem. In other words,

Instead of

    Hubbard sites
    if( sitesfile exists)
       {   
       readFromFile(sitesfile,sites);
       }
    else
      {
       sites = Hubbard(N_sites);
       auto sites1 = sites;
       auto sites2 = sites;
       }


I used

    Hubbard sites, sites1, sites2;
    if( sitesfile exists)
      {
      readfromfile(sitesfile,sites);
      }
     else
      {
       sites = Hubbard(N);
       sites1 = sites;
       sites2 = sites;
      }

---

Yes, I hadn't noticed before that you had two different "sites" variables sites1 and sites2. I would discourage doing this because the main idea of a site set is to ensure that various MPS all have the same site indices. So the intended use of a site set is that you have a single site set in your program that you use to initialize all of your MPS and your InitState objects. In your code above things work ok because sites1 and sites2 are just copies of the same original site set. But ideally you would just use a single site set throughout.

---

Ok thanks for the reminder, I will try it that way to stick with the intended use of sites variables.

---

Hi Miles,

Thanks for reminding me. When I normalized psi, I got 122...2 without "psi.position(j)" and 11...1 with "psi.position(j)".