An non-interacting fermion with nonzero flux

Question

Hi Miles,

I encountered a problem recently when I tried to implement the hopping part of the extended Hubbard model with nonzero flux. The following is a part of my code:
if(phi==0)
{
ampo += -tpara,"Cdagup",b,"Cup",b+2;
ampo += -tpara,"Cdagup",b+2,"Cup",b;
ampo += -tpara,"Cdagdn",b,"Cdn",b+2;
ampo += -tpara,"Cdagdn",b+2,"Cdn",b;
}
else
{
ampo += (-(tpara)std::polar(1.,-phi/4.)),"Cdagup",b,"Cup",b+2;
ampo += (-(tpara)std::polar(1.,+phi/4.)),"Cdagup",b+2,"Cup",b;
ampo += (-(tpara)std::polar(1.,-phi/4.)),"Cdagdn",b,"Cdn",b+2;
ampo += (-(tpara)std::polar(1.,+phi/4.)),"Cdagdn",b+2,"Cdn",b;

}

The variable phi indicates the flux. The strangest thing is that this Hamiltonian has the same energy as the result obtained by the exact diagonalization with 0.5*phi flux, but it should not. It should give the result with phi/4 flux. I will be very thankful if you can help me to find out where the problem is in my code.

Comments

Also I used the exactly the same code but with ITensor 1 and found different results which are the same as the exact diagonization. It is quite confusing for me.

Answer 1

Hi, so I tested a model very similar to the wrote you wrote above and AutoMPO seems to be making it correctly.

To see the test I wrote, do the following:
1. Do 'git pull' to get the latest version of ITensor
2. Look in the unittest folder, at the file unittest/autompo_test.cc
3. At the end of the file, you'll see a test called "Hubbard, Complex Hopping"

You can use similar code to test your own Hamiltonians. The key technique is to use the InitState class to make various product state wavefunctions, then use the functions overlap or overlapC to compute matrix elements of your Hamiltonian. For instance, you can use states containing a single particle to test hopping terms.

I would not fully trust the version 1 AutoMPO for fermions at the moment. I have a new version of it in development (on the v1 branch) that I should push soon, but for now please use version 2 which I have been keeping much more up-to-date.

Finally, there is a suspicious looking thing about the Hamiltonian you wrote: it only connects odd sites to odd sites and even sites to even sites. That is, it divides the system into two sublattices which don't couple to one another. Could this perhaps be responsible for the incorrect results you saw?

Miles

Comments

Hi, I am thankful for your reply. However my problem remains even I have updated my code. Actually my Hamiltonian is a two-leg hubbard model. There are two types of hopping. One is along the parallel direction. Another hopping is between the leg. The parallel hopping has extra phase to represent the flux and the hopping between leg is real. Therefore my Hamiltonian has both real and complex elements. The following is my full kinetic part:

for(int b=1; b<=N-3; b+=2)
{
    if(phi==0)
    {
      ampo += -tpara,"Cdagup",b,"Cup",b+2;
      ampo += -tpara,"Cdagup",b+2,"Cup",b;
      ampo += -tpara,"Cdagdn",b,"Cdn",b+2;
      ampo += -tpara,"Cdagdn",b+2,"Cdn",b;
    }
    else
    {
      ampo += (-(tpara)*std::polar(1.,-phi/2.)),"Cdagup",b,"Cup",b+2;
      ampo += (-(tpara)*std::polar(1.,+phi/2.)),"Cdagup",b+2,"Cup",b;
      ampo += (-(tpara)*std::polar(1.,-phi/2.)),"Cdagdn",b,"Cdn",b+2;
      ampo += (-(tpara)*std::polar(1.,+phi/2.)),"Cdagdn",b+2,"Cdn",b;
    }
    ampo += -tperp,"Cdagup",b,"Cup",b+1;
    ampo += -tperp,"Cdagup",b+1,"Cup",b;
    ampo += -tperp,"Cdagdn",b,"Cdn",b+1;
    ampo += -tperp,"Cdagdn",b+1,"Cdn",b;
  }
  ampo += -tperp,"Cdagup",N-1,"Cup",N;
  ampo += -tperp,"Cdagup",N,"Cup",N-1;
  ampo += -tperp,"Cdagdn",N-1,"Cdn",N;
  ampo += -tperp,"Cdagdn",N,"Cdn",N-1;
  for(int b=2; b<=N-2;b+=2)
  {
    if(phi==0)
    {
      ampo += -tpara,"Cdagup",b,"Cup",b+2;
      ampo += -tpara,"Cdagup",b+2,"Cup",b;
      ampo += -tpara,"Cdagdn",b,"Cdn",b+2;
      ampo += -tpara,"Cdagdn",b+2,"Cdn",b;
    }
    else
    {
      ampo += (-(tpara)*std::polar(1.,+phi/2.)),"Cdagup",b,"Cup",b+2;
      ampo += (-(tpara)*std::polar(1.,-phi/2.)),"Cdagup",b+2,"Cup",b;
      ampo += (-(tpara)*std::polar(1.,+phi/2.)),"Cdagdn",b,"Cdn",b+2;
      ampo += (-(tpara)*std::polar(1.,-phi/2.)),"Cdagdn",b+2,"Cdn",b;
    }
  }
  IQMPO H = IQMPO(ampo);

I set N=12 and number of particles is 6. The exact solution with phi=\pi is -11.4401 but I can
not reproduce with my code. Actually my code gives me the solution the same as the solution with zero flux. When I turn off the flux, the code behave normally.  Could you help me to solve my problem?

---

Based on what you're reporting, it's hard to know whether the issue is with H or with DMRG (or possibly the initial state wavefunction). It could be any of these things: for instance, you may need to start with a complex wavefunction, or turn on noise during the initial sweeps.

One thing it would be a very good idea to do, though, is to check the Hamiltonian matrix elements. You can use code like this:

        auto state1 = InitState(sites,"Emp");
        state1.set(j,"Up");
        auto state2 = InitState(sites,"Emp");
        state2.set(j+1,"Up");
        auto psi1 = IQMPS(state1);
        auto psi2 = IQMPS(state2);
        Print(overlap(psi2,H,psi1) - z);

where z is the expected matrix element of H that hops a fermion from site j to j+1. If the result is zero you can be sure the Hamiltonian is defined correctly. Then we'll know that is has something to do with the DMRG and not AutoMPO.

---

Hi,
The following is my test part
for(int i=1;i<=N-3;i+=2)
  {
        auto state1 = InitState(sites,"Emp");
        state1.set(i,"Up");
        auto state2 = InitState(sites,"Emp");
        state2.set(i+2,"Up");
        auto psi1 = IQMPS(state1);
        auto psi2 = IQMPS(state2);
        Print(overlapC(psi2,H,psi1) - (-(tpara)*std::polar(1.,+phi/2.)));
  }
  for(int i=2;i<=N-2;i+=2)
  {
        auto state1 = InitState(sites,"Emp");
        state1.set(i,"Up");
        auto state2 = InitState(sites,"Emp");
        state2.set(i+2,"Up");
        auto psi1 = IQMPS(state1);
        auto psi2 = IQMPS(state2);
        Print(overlapC(psi2,H,psi1) - (-(tpara)*std::polar(1.,-phi/2.)));
  }
The first loop is correct but the second is wrong when phi=\pi. However when I turn on an arbitrary flux. Both loops are wrong. I suspect the problem probably in the conversion of AutoMPO to IQMPO.

---

Ok, thanks a lot for checking. Since it seems to be an AutoMPO issue, I will run your test code, see if I can reproduce, and if so fix it.

---

Ok good news / bad news. The good news is that I could get the AutoMPO conversion to work properly by replacing the line H = IQMPO(ampo) with this line:
auto H = toMPO<IQTensor>(ampo,{"Exact",true});

This calls a different backend engine for generating the IQMPO which is less complicated and gives the correct result when I did a matrix element test.

The bad news (for me) is that you are definitely right about there being a bug in the current default AutoMPO conversion engine. Thanks for finding it and your patience with checking the Hamiltonian and everything.

Please give that new approach a try (the one with auto H = toMPO<IQTensor>(ampo,{"Exact",true}); ) and let me know if it still doesn't give you the right energy.

---

Hi Miles,
Thank you for your suggestion. This approach gives me the correct matrix element and energy for my two-leg hubbard model with and without flux. But I am now wondering why you don't use this engine as default at beginning. Anyway, thank you again to solve my problem.

---

The "Exact" engine that is currently working is the old engine, which is simpler and makes no approximations, but has the limitation of only allowing at most two operators in each term.

The new engine lets you have more than 2 operators in a term and can also deal with long-range interactions much better. But it's a lot more complicated, thus prone to be more buggy at the moment.

---

Ok, so I found the bug pretty quickly! There is a part where we do an SVD then apply V^\dagger to certain matrices to construct the MPO/IQMPO. The conjugate part of the \dagger operation was missing in some key places.

I just pushed the fixed version which also includes a thorough unit test of your ladder Hamiltonian to make it impossible for this bug to occurr again in the future (without the test failing I mean).

Thanks for your patience and for finding this bug. I'll credit you when I write the changelog release notes for the next version.