How to realize Hubbard-Holstein model

Question

I want to realize Hubbard-Holstein model which is mixed bose-fermi model. I watch the question:enter link description here

the QN file work well. However, Hubbard-holstein hamiltonian like below(Cdagup,Cdagdn,Cup,Cdn is spinfull fermi and Ab,Adagb is bose). I want to initial fermi state and bose state. When I set the initial state , the program error: setting IQtensor element non-zero would violate its symmetry. I think the problem is

    ampo += g,"Cdagup",j,"Ab",j+1;

which are not conserve bose number. I want to know how to correct code,Thank a lot.

int N=80;
auto U=3.6;
auto J=1.0;
auto g=3.0;
auto K=5.0;
auto sites=BF(N);
auto ampo=AutoMPO(sites);
for(int j=1;j<=N-2;j+=2)
    {
    ampo += -J,"Cdagup",j+2,"Cup",j;
    ampo += -J,"Cdagdn",j+2,"Cdn",j;
    ampo += -J,"Cdagup",j,"Cup",j+2;
    ampo += -J,"Cdagdn",j,"Cdn",j+2;
}

for(int j=1;j<=N;j=j+2)
{
ampo += U,"Nup",j,"Ndn",j;
ampo += g,"Cdagup",j,"Ab",j+1;
ampo += g,"Cup",j,"Adagb",j+1;
ampo += g,"Cdagdn",j,"Ab",j+1;
ampo += g,"Cdn",j,"Adagb",j+1;
ampo += K,"Adagb",j+1,"Ab",j+1;
}

}

auto H = IQMPO(ampo);
auto state0 = InitState(sites);
for(int i = 1; i <= N; i=i+1) 
{
    if (i%4==1)  
    {
        state0.set(i,"Up");
    }
    else if (i%4==3)   
    {
        state0.set(i,"Dn");
    }
    else 
    {
        state0.set(i,"b1");
    }
}

auto psi = IQMPS(state0);

auto sweeps=Sweeps(5);
sweeps.maxm() = 20,50,100,200,500;
sweeps.cutoff() = 1E-12;


auto energy = dmrg(psi,H,sweeps,{"Quiet=",true});

Comments

Hi Weijie,
Could you please post the exact (mathematical) description of the Hamiltonian you are trying to make? From the literature, I believe in the Holstein model the "Ab" (boson annihilation) operator couples to the fermion density operator, not to the fermion creation/annihilation operator as in the example AutoMPO term you showed.

Best regards,
Miles

---

Dear miles,

I'm sorry I write it wrong. The hamiltonian
ampo += g,"Cdagup",j,"Ab",j+1;
ampo += g,"Cup",j,"Adagb",j+1;
ampo += g,"Cdagdn",j,"Ab",j+1;
ampo += g,"Cdn",j,"Adagb",j+1;
should be corrected as:
ampo += g,"Nup",j,"Ab",j+1;
ampo += g,"Nup",j,"Adagb",j+1;
ampo += g,"Ndn",j,"Ab",j+1;
ampo += g,"Ndn",j,"Adagb",j+1;

The hamiltonian which I want to realize is
H=-t\sum_{i\sigma}(c_{i+1,\sigma}^{\dag}c_{i,\sigma}+h.c.)+U\sum_{i}n_{i\uparrow}n_{i\downarrow}+g\sum_{i,\sigma}n_{i,\sigma}(b_{i}+b_{i}^{\dag})+K\sum_{i}b_{i}^{\dag}b_{i}. b^{\dag} and b represent boson and the other is fermion.

Thanks for your quick reply,
Weijie Huang

Answer 1

Hi Weijie,
Ok thanks, that looks more like what I expected. So to implement this model, you will have to make a custom site set which is not one of the standard ones provided. But it should be easy to make.

Steps to do:

1. upgrade to ITensor version 3 if you haven't already:
   http://itensor.org/docs.cgi?vers=cppv3&page=upgrade2to3

2. define your new site set "Holstein" as follows (outside of your main function):

   using Holstein = MixedSiteSet<ElectronSite,BosonSite>;

3. now if you make a site set object like this:

   auto sites = Holstein(N,{"MaxOcc=",3});

then the odd numbered sites will be electron sites and the even numbered sites will be boson sites. The "MaxOcc" named argument is what sets the maximum occupancy of the boson sites. Of course in reality one would like this to be very large, but practically one has to limit this to a finite value; you should check that your results are converged in this value and that you have taken it large enough.

1. when you make your AutoMPO, be sure to use only operators defined for the Electron site set (http://itensor.org/docs.cgi?vers=cppv3&page=classes/electron) on odd-numbered sites and the boson site set (http://itensor.org/docs.cgi?vers=cppv3&page=classes/boson) for even numbered sites

2. once your Hamiltonian gets made correctly (possibly after further discussion here), make sure to test it carefully on small systems and check, if possible, with results obtained using other methods, or at least in limits where you have a good understanding of what the results ought to be

---

Following up on the comments below, here is a complete working sample code (to be used with the very latest version of the v3 branch, dated May 29):

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

using Holstein = MixedSiteSet<ElectronSite,BosonSite>;

int
main()
    {
    auto N = 20;

    auto U = 1.0;
    auto J = 1.0;
    auto g = 3.0;
    auto K = 5.0;

    auto sites = Holstein(N,{"ConserveNf=",true,
                             "ConserveNb=",false,
                             "MaxOcc=",3});

    auto ampo = AutoMPO(sites);
    for(int j=1;j  <= N-2; j+=2)
        {
        ampo += -J,"Cdagup",j+2,"Cup",j;
        ampo += -J,"Cdagup",j,"Cup",j+2;
        ampo += -J,"Cdagdn",j+2,"Cdn",j;
        ampo += -J,"Cdagdn",j,"Cdn",j+2;
    }
    for(int j=1;j < N; j += 2)
        {
        ampo += U,"Nup",j,"Ndn",j;
        ampo += g,"Ntot",j,"A",j+1;
        ampo += g,"Ntot",j,"Adag",j+1;
        }
    for(int j = 2; j <= N; j += 2)
        {
        ampo += K,"N",j;
        }
    auto H = toMPO(ampo);

    auto sweeps = Sweeps(20);
    //very important to use noise for this model
    sweeps.noise() = 1E-6,1E-6,1E-8,1E-12;
    sweeps.maxdim() = 10,20,100,100,200;
    sweeps.cutoff() = 1E-10;

    auto state = InitState(sites);
    for(auto j : range1(N))
        {
        if(j%4==1) state.set(j,"Up");
        else if(j%4==3) state.set(j,"Dn");
        else state.set(j,"0");
        }
    auto psi0 = MPS(state);

    auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet=",true});

    printfln("Ground State Energy = %.12f",energy);

    auto checkNtot = 0.0;

    println("Up Electron Density:");
    for(int j = 1; j <= N; j += 2)
        {
        psi.position(j);
        auto upd = elt(dag(prime(psi(j),"Site"))*op(sites,"Nup",j)*psi(j));
        checkNtot += upd;
        printfln("%d %.12f",j,upd);
        }
    println();

    println("Dn Electron Density:");
    for(int j = 1; j <= N; j += 2)
        {
        psi.position(j);
        auto dnd = elt(dag(prime(psi(j),"Site"))*op(sites,"Ndn",j)*psi(j));
        checkNtot += dnd;
        printfln("%d %.12f",j,dnd);
        }
    println();

    Print(checkNtot);

    println("Boson Density:");
    for(int j = 2; j <= N; j += 2)
        {
        psi.position(j);
        auto bd= elt(dag(prime(psi(j),"Site"))*op(sites,"N",j)*psi(j));
        printfln("%d %.12f",j,bd);
        }
    println();

    return 0;
    }

Comments

Dear miles,

I take some time to read the new change of version 3. It seems to be more convenient. According to your reply , I change my codes as:

#include "itensor/all.h"

using namespace itensor;
using Holstein = MixedSiteSet<ElectronSite,BosonSite>;

int
main()
    {
    int N = 20;
    auto U=3.6;
    auto J=1.0;
    auto g=3.0;
    auto K=5.0;
    auto sites = Holstein(N,{"ConserveNf",true,"MaxOcc=",3});

    auto ampo = AutoMPO(sites);
    for(int j=1;j<=N-2;j+=2)
        {
        ampo += -J,"Cdagup",j+2,"Cup",j;
        ampo += -J,"Cdagdn",j+2,"Cdn",j;
        ampo += -J,"Cdagup",j,"Cup",j+2;
        ampo += -J,"Cdagdn",j,"Cdn",j+2;
    }
   for(int j=1;j<=N;j=j+2)
   {
        ampo += U,"Nup",j,"Ndn",j;
        ampo += g,"Ntot",j,"A",j+1;
        ampo += g,"Ntot",j,"Adag",j+1;
        ampo += K,"Adag",j+1,"A",j+1;
        }
    auto H = toMPO(ampo);

    auto sweeps = Sweeps(5);
    sweeps.maxdim() = 10,20,100,100,200;
    sweeps.cutoff() = 1E-10;

    auto psi0 = randomMPS(sites);

    auto [energy,psi] = dmrg(H,psi0,sweeps);

    println("Ground State Energy = ",energy);

    return 0;
    }

the compile is successful,but finally there is a error: Index does not contain given QN block. I don't know what is problem.

Best regards,
Weijie Huang

---

Dear miles,

Maybe I add the button to select the answer, so you don't answer me. However, my problem is not solved. I hope get the answer.

Best regards,
Weijie Huang

---

Hi Weijie, thanks for your comment. For some reason I didn't get an email about your previous comment, so I'm glad you posted a new one to remind me to respond to you.

It turns out the reason for you bug was kind of subtle, but we made the necessary changes to ITensor to make the Hubbard-Holstein model work correctly. One has to pass the named args {"ConserveQNs=",true,"ConserveNb=",false} to get that model to work, since it does not conserve boson number by definition. Also you'll note in the sample code above there is a MaxOcc named argument which you should also adjust to make sure it is not affecting the results you get.

Please see my edited answer above with a complete sample code for simulating the model. Also please note the importance of including the "noise" setting of the sweeps which turns on the "DMRG noise term" which is very crucial for converging complicated models of this type with multiple quantum numbers and complicated interactions.

Finally note that the way the lattice is defined it has a left-right asymmetry since the first site is an electron site but the last site is a boson site. This is ok, and not a problem, it's just that you may want to take this into account when interpreting your results, or define the lattice or Hamiltonian slightly different. But of course it shouldn't matter in the limit of a very long chain.

Finally please make sure to do a large number of sweeps for a model like this to ensure that it is converged. Experiment with different numbers of sweeps to make sure.

Best regards,
Miles

---

Dear miles

I only substitude my code :auto sites = Holstein(N,{"ConserveNf",true,"ConserveNb=",false}); for       auto sites = Holstein(N,{"ConserveNf",true,"MaxOcc=",3});
However, the Hubbard-Holstein model don't work correctly. The error is: Index does not contain given QN block.The error seems same as before.

Best regards,
Weijie Huang

---

Dear miles

Thanks for your help, when I update the latest version. The calculation is correct.

Best regards
Weijie Huang