Questions about spin and charge gap on Holstein-Hubbard model

Question

I'm following an article: arXiv:0808.1675. When I calculate the spin and charge gap of Fig.2, I can't get the correct answer,my codes is:

#include "itensor/all.h"

using namespace itensor;

int
main()
{
using Holstein =MixedSiteSet<ElectronSite,BosonSite>;
auto N=60;
auto U=1.0;
auto t=1.0;
auto eplison=1.25;
auto omega=5.0;
auto sites=Holstein(N,{"ConserveNf",true,"ConserveNb",false,"MaxOcc=",4});

auto ampo =AutoMPO(sites);
for(int j=1;j<=N-2;j=j+2)
{
    ampo += -t,"Cdagup",j+2,"Cup",j;
    ampo += -t,"Cdagdn",j+2,"Cdn",j;
    ampo += -t,"Cdagup",j,"Cup",j+2;
    ampo += -t,"Cdagdn",j,"Cdn",j+2;
}
for(int j=1;j<=N;j=j+2)
{
    ampo += U,"Nup",j,"Ndn",j;
    ampo += -sqrt(eplison*omega),"Nup",j,"A",j+1;
    ampo += -sqrt(eplison*omega),"Nup",j,"Adag",j+1;
    ampo += -sqrt(eplison*omega),"Ndn",j,"A",j+1;
    ampo += -sqrt(eplison*omega),"Ndn",j,"Adag",j+1;
    ampo += omega,"Adag",j+1,"A",j+1;
}
auto H=toMPO(ampo);

auto sweeps =Sweeps(20);
sweeps.noise()=1E-6,1E-6,1E-8,1E-12;
sweeps.maxdim()=10,20,100,100,200,400,800;
sweeps.cutoff()=1E-10;

auto state =InitState(sites);
for(auto j:range1(N))
{
    if(j%4==1)state.set(j,"Dn");
    else if(j%4==3)state.set(j,"Up");
    else state.set(j,"0");
}
auto psi0=randomMPS(state);
auto [energy,psi]=dmrg(H,psi0,sweeps,{"Quiet=",true});


auto state1=InitState(sites);
for(auto j:range1(N))
{
    if(j==1)state1.set(j,"UpDn");
    else if(j%4==1&&j!=1)state1.set(j,"Dn");
    else if(j%4==3)state1.set(j,"Up");
    else state1.set(j,"0");
}
auto psi1=randomMPS(state1);
auto [energyadd,psi2]=dmrg(H,psi1,sweeps,{"Quiet=",true});


auto state2=InitState(sites);
for(auto j:range1(N))
{
    if(j==3)state2.set(j,"Emp");
    else if(j%4==1)state2.set(j,"Dn");
    else if(j%4==3&&j!=3)state2.set(j,"Up");
    else state2.set(j,"0");
}
auto psi3=randomMPS(state2);
auto [energyminus,psi4]=dmrg(H,psi3,sweeps,{"Quiet=",true});

auto state3=InitState(sites);
for(auto j:range1(N))
{
    if(j==1)state3.set(j,"Up");
    else if(j%4==1&&j!=1)state3.set(j,"Dn");
    else if(j%4==3)state3.set(j,"Up");
    else state3.set(j,"0");
}
auto psi5=randomMPS(state3);
auto [energyspin,psi6]=dmrg(H,psi5,sweeps,{"Quiet=",true});
printfln("spin gap=",energyspin-energy);
printfln("charge excitation=", energyadd+energyminus-2*energy);
    return 0;
   }

When I calculate N=60 in the code(actually the system is 30) , lambda=0.625, boson is 4. the spin gap is 0.351311, the charge gap is 0.442756. They are not correct. I try the other ways, but it doesn't work. Whether I ignore something?

Comments

Hi Weijie,
At a glance your code looks correct, but it’s hard to know what is happening without running it myself, and these kind of questions about details of a particular run are difficult to answer. Did you try running your code on a very small system and comparing to exact results to check it for correctness? That’s always a good idea and something I often do with a new code.

Then if you have tested your code is logically correct, when using it to get excited states for larger systems I’ve found a number of different details can matter quite a lot, including:
- the number of sweeps you do (it needs to be rather large, like 50 or 100 often for excited states)
- how well converged the ground state is, against which you are orthogonalizing the higher-energy states
- the choice of initial state(s) used to initalize the DMRG calculations (this is important to experiment with)

There is always the chance that our excited-state DMRG code will not return the excited states in a strictly sequential manner. For example it could return the second excited state when you provide the ground state, then when you provide those two states it might next give you the first excited state. It is a tricky and more complicated optimization than just finding the ground state.

Also please pay attention to all of the quantum numbers you set (through your choice of initial state). Are you sure the spin gap is to a state with the same quantum numbers as the excited state you calculated?

Hope that helps -

Miles

---

Thanks for your patient answer.

For the excited state. Actually I only calculate ground-state. the charge and spin gaps are determined from:
\Delta_{c}=E_{0}^{+}(1/2)+E_{0}^{-}(-1/2)-2E_{0}(0)
\Delta_{s}=E_{0}(1)-E_{0}(0)
where \E_{0}^{\pm}(S^{z}) is the ground-state energy at(away from) half-filling with N_{e}=N(N_{e}=N\pm 1)

For the initial state, I have tried change the randomMPS to MPS to find answer. However, it wasn't right.

I want to know if I should update version? I don't know if this influence the answer.
Best wishes

---

Hi Weijie,
Thanks for clarifying that - I see now that's what the code is doing.

So there are a lot of possibilities and it's hard to answer without hearing more information from you. What are the correct values of the gaps that you expected, and how different are they from the numbers you got?

Also when you run the DMRG calculations and watch the energy at each sweep, is it still going down in the last 4 or 5 sweeps or is it converged? Does increasing the maximum allowed bond dimension or lowering the DMRG cutoff help?

Also, I notice you are working on somewhat small system sizes of 30 unit cells. What happens when you do 40 instead? Or 50? Do the gaps get a lot closer to the values you expect? Or are those values themselves finite-size values?

Finally, no I don't think updating ITensor is the most likely reason, as our DMRG code does not usually change often and it gives pretty stable results across versions. But at the same time, we encourage users to upgrade to the latest tagged version so that they can benefit from bug fixes we put in.

Best regards,
Miles

---

I have tried system size  30, spin gap is 0.265792( the value I expect is around 0.33), in the last 4 sweeps, energy are -80.136093617863, -80.136094691090, -80.136094766538, -80.136094763304
        system size is 40,spin gap is 0.203344(the value I expect is around 0.3),in the last 4 sweeps, energy are -107.183336662006, -107.183339581477, -107.183339990703, -107.183340011832.
        system size is 50, spin gap is 0.167828(the value I expect is around 0.25), in the last 4 sweeps,energy are -134.206961626752, -134.207171511814, -134.207261444963, -134.207280146347.
        It seems that the gaps don't get closer when the system is larger.
Best regards,
miles

---

At a glance, it seems that this system has multiple, degenerate ground states. Have you taken that into account when choosing which energies to subtract to define the gap?

Miles

---

Hi Miles,
        Actually the article doesn't tell me  how to choose, if there are many degenerate ground states, I have to read more papers. However, I want to know why do you think it has degenerate states? And whether there are ways to find degenerate states?
Best regards,
Weijie Huang

---

Hi Weijie, I see - sorry. I had misunderstood the numbers you posted as being the energies of different states you found, but now I see they are just the energies in the last 4 sweeps so of course they are very close together.

But now I’m confused about how the information you showed indicates anything about the gaps or that they are getting closer together. To see how the gap is behaving, I would need to see the energies of two different states. Are the last-four-sweep energies you are showing for the ground state or the first excited state you’re using to define the gap?

Miles

---

Hi, Miles
Actually I only calculate ground-state. I want to calculate spin gaps which are determined from:
\Delta_{s}=E_{0}(1)-E_{0}(0)
where E_{0}(1) is the ground state at half-filled with sz=1, and E_{0}(0) is the ground state at half-filled with sz=0
Weijie Huang

---

Ok I see. I recall that now. But in your post above you only listed one set of energies for a given system size, so that confused me. Of course you'll be looking at two different energies for each system size, as you know, both of which are ground states in their own sector.

So if the difference between these two energies does not get close to the expected value for the gap, it could be for a number of different reasons and it's hard to know which without doing more tests and comparisons. Here are some possible reasons, in order of most likely to least likely:
1. most likely one or both of your calculations are stuck in a metastable or excited state rather than finding the true ground state for that quantum number sector. This can happen for challenging Hamiltonians with DMRG (not specific to ITensor, but to the DMRG algorithm). The way to fix this is to use a better guess for an initial state, though a randomMPS with a large bond dimension could be a very good guess. Also you can turn on the noise term if you aren't already.

2. it could be that there is a large finite-size effect in your system or that you are otherwise not comparing exactly the same quantity to what you'd expect from a different calculation, like if it uses different boundary conditions etc.

3. it could be that there is a bug in your code, such as one of the MPS not being set to the proper quantum number sector you expect it to, or an error in the definition of your Hamiltonian

So it's hard to say from just looking at the code and some of the results. I'd recommend testing each of these above 3 things more in your code by printing out more details such as the totalQN of each MPS, or taking matrix elements of product states with your Hamiltonian MPO to check that it's correct, and inputting other initial states into DMRG to see if you get a lower energy etc.

Hope that helps -

Miles

---

Hi, miles
      Thanks for your answer,  I'll try these suggestions. If there are some more information, I'll post it.
Best ragard,
Weijie Huang