Excited states in the Hubbard model

Question

Hello all,

I have been trying to get excited states of the Hubbard model as an approximation to the continuum.

I have essentially adapted the example provided here:

https://www.itensor.org/docs.cgi?vers=cppv3&page=formulas/excited_dmrg

Here is the code:

#include "itensor/all.h"

using namespace itensor;

int main()
{
auto N = 200;//system size
auto Npart = 5; //particle number
auto L = 10.;
float a = L/N;
float t1 = 0.5/(aa);
auto g = 10.;
float U = g/a;
auto Ns = 10; //number of excited states
float s = 1/(aa);

println("number of sites ",N);
println("system length ",L);
println("lattice spacing ",a);
println("hopping ",t1);

auto sites = Electron(N,{"ConserveQNs=", true});
auto ampo = AutoMPO(sites);

for(int i = 1; i <= N; ++i) 
    {
    ampo += U,"Nupdn",i;
    }
for(int i = 1; i <= N; ++i) 
    {
    ampo += s,"Nup",i;
    ampo += s,"Ndn",i;
    //ampo += V*(i),"Nup",i;
    //ampo += V,"Ndn",i;
    } 
for(int b = 1; b < N; ++b)
    {
    ampo += -t1,"Cdagup",b,"Cup",b+1;
    ampo += -t1,"Cdagup",b+1,"Cup",b;
    ampo += -t1,"Cdagdn",b,"Cdn",b+1;
    ampo += -t1,"Cdagdn",b+1,"Cdn",b;
    }
auto H = toMPO(ampo);


auto state = InitState(sites);
int p = Npart;
for(int i = N; i >= 1; --i) 
    {
    if(p > i)
        {
        println("Doubly occupying site ",i);
        state.set(i,"UpDn");
        p -= 2;
        }
    else
    if(p > 0)
        {
        println("Singly occupying site ",i);
        state.set(i,(i%2==1 ? "Up" : "Dn"));
        p -= 1;
        }
    else
        {
        state.set(i,"Emp");
        }
    }

auto psiA = MPS(state);

auto sweeps = Sweeps(25);
sweeps.maxdim() = 10E4;
sweeps.cutoff() = 1E-10;
sweeps.niter() = 4;
sweeps.noise() = 1E-7,1E-8,0.0;
println(sweeps);

auto [en0,psi0] = dmrg(H,psiA,sweeps,"Quiet");

auto wfs = std::vector<MPS>(0);
wfs.push_back(psi0);
Vector energies(Ns+1);
energies(0) = en0;

for(int j = 1; j <= Ns; ++j)
    {
    auto [en,psi] = dmrg(H,wfs,psi0,sweeps,{"Quiet=",true,"Weight=",100.0});
    energies(j) = en;
    wfs.push_back(psi);
    printfln("\n%fth Excited State Energy = %f",j,en);
    printfln("\nOverlap <psi0|psi%f> = %E",j,inner(psi0,psi));
    }

std::ofstream file("hub_spectrum.dat");

for(int j = 0; j <= Ns; ++j)
    {
    printfln(file,"%d  %f",j,energies(j));
    }

file.close();

return 0;
}

I have a few issues when running this:

• I If I try to use psiA (the Neel state) as the initial guess for the excited state search I get the error

div(T1)=QN({"Sz",-2}) must equal div(T2)=QN() when adding T1+T2

which has also been reported here:

http://itensor.org/support/1581/error-when-calculating-excited-states-conserved-quantities?show=1581#q1581

• If I use psi0 (the ground state found by dmrg) as the initial guess, I get a bunch of energies, but they are often out of order (not necessarily increasing in value). In fact, the first excited state seem nearly degenerate with the ground state, which should not be the case. Notice that I am requiring that the new state be orthogonal to the previous ones.

So my questions are:

1. How can I optimise the code to get the correct excited state energies, in the right order?
2. Is there a particular choice of initial state guess that would solve these issues (and why does that neel state fails for the excited state search)?
3. Making the "Weight" parameter larger does not necessarily improve the result. Why is that?

Best regards,

Rafael

Answer 1

Hi Rafael,
This is similar to a bug some other users have reported. I just patched ITensor to hopefully fix it. Could you please "git pull" the latest version of ITensor on the v3 branch, fully recompile the library, then check if your code now works as expected? Thanks!

Miles

Comments

Dear Miles,

Thanks for the reply. I have done as you suggested and indeed that seems to solve the problem! I will keep running a few tests and I'll get back to you if there's anything else.

Do you have any ideas on the other points (appropriate choice of initial states and weight parameter) or this should not have a big influence on the outcome?

Best wishes,

Rafael

---

Dear Rafael,
Those are good questions, but they are also hard to answer in some fully complete way, and in a way are kind of empirical and system-dependent.

For initial states, it really depends on the physics, and the nature of the states you are trying to find. It's hard to be more specific than that. Often it's good enough to just take a product state with the right total QN and which is not too pathological; e.g. for particles making sure the particles are well spread out over the system. But if you suspect the initial state is causing problems and DMRG is getting stuck, then various things can be tried such as doing a small amount of imaginary time evolution before switching to DMRG, using the ground state to build a guess for the excited state(s); using randomized states; etc.

For the weight value, it's also not totally clear. It's one of those things where the best answer in practice is to just try various values of it and check that the results don't depend strongly on which value you choose. The weight parameter is the coefficient of a term added to the Hamiltonian which penalizes any overlap between the state being optimized and the previous states passed to DMRG. So if it's too small then the overlap may not go to a small enough value at the end; too large and DMRG may struggle to converge or find the wrong state etc. Really it's just important to explore different ranges of this parameter for your system, and it depends on the sizes of the gaps between states and other things probably.

Hope that gives you some good ideas and starting points -

Miles

---

Hi Miles,

I think these are good guidelines, thanks for the reply.

Best regards,

Rafael