+1 vote
asked by (360 points)

Hello everyone!

I have been trying to use ITensor to calculate the properties of a Hubbard Hamiltonian with 3 internal components.

I think the questions is related in a way to this one:

http://itensor.org/support/25/multiband-hubbard-model?show=25#q25

although I cannot see how to construct such a model by simply extending the Hubbard code with more sites in an unit cell.

How can we construct the set of operators required to approach this problem? Do we have examples of that?

Best,

Rafael

1 Answer

+2 votes
answered by (70.1k points)

Hi Rafael,
I think your idea of using more than one site in a unit cell is the right thing to do for creating multi-component / multi-band models of interacting electrons. Basically let's say you have 3 "flavors" or "labels" of electrons. Then sites 1,4,7,10 etc. are all considered (by you) to be flavor 1; sites 2,5,8, etc. flavor 2; and 3,6,9, etc. flavor 3.

You then just define the Hamiltonian parameters as from on paper, so like if t1 is the hopping between sites which are both flavor 1, you just have t1 go between every third site, being sites 1,4,7, (i.e. the flavor 1 sites). That is it's a 3rd neighbor hopping. Similar for a t2.

In contrast, a local term like a Hunds interaction would happen inside the unit cell, so like between just sites 2 and 3, then sites 5 and 6, etc. repeating that way.

If that's still not too clear, perhaps you could provide the specific Hamiltonian or a link to it and we could discuss more?

Best,
Miles

commented by (360 points)
Hi Miles,

Thanks for the reply! Sorry for not writing the detailed Hamiltonian, but I was thinking more in terms of a general three-component model.

But that is indeed a very clever way to do it, thanks for the details. I am assuming that you are considering a single-component or spinless fermion model as your starting point there, right?

I will try to implement that and I'll come back if there are more questions.  

Best wishes,

Rafael
commented by (70.1k points)
Hi Rafael,
Sure thing. To expand a bit more, if the Hamiltonian you are considering is a multi-band model of electrons, then having more than one site (degree of freedom) in the unit cell is essentially the defining case of how multiple bands arise in real systems. So far from being a trick, it's actually the closest to reality in some sense.

Regarding whether the fermions are spinless or have spin, I can't answer that because I don't know what model you are studying. Do the fermions in your model have spin? That's something I should be asking you, instead of you asking me :^)

Best,
Miles
commented by (360 points)
Hi Miles,

Thanks again. I think I was not clear before: I am considering a three component fermionic system such as the ones that are realised nowadays with cold atoms:

https://arxiv.org/pdf/1502.02495.pdf

Eq. 1 shows the type of Hamiltonian that I am interested in. So it is indeed a fermionic system with "spin" in the sense that these are internal states.

Now, what I mean by having the spinless fermion system as a starting point: since you are suggesting to represent the internal states as lattice sites, then I am assuming that the set of operators is created with the Fermion class, that is, no spin states. That is what makes sense to me.

Now, here is another question: is there a way to guarantee the conservation of particle number in the different components with such an approach? Clearly the total fermion number could be conserved, but I see some situations (in the case of asymmetrical interactions between different components, for instance) where this might not be the case.

Best wishes,

Rafael
commented by (70.1k points)
Hi Rafael, I see - that's very helpful to have the paper to look at.

To answer your questions:

(1) yes, I'd recommend using Fermion sites. (Technically you could use Electron sites but it might be kind of awkward to map the spin labels to component labels even though it could be done in some cases.) But in a certain case, I'd actually recommend you make your own modified version of Fermion sites to use - see the next part of the answer below

(2) if you are definitely wanting to conserve the different components (different \alpha labeled fermions), then what you'll need to do is to define distinct quantum number names. So instead of "Nf" which is our default convention for "number of fermions" you could use new names like "Nfa", "Nfb", "Nfc" or "Nf1", "Nf2", "Nf3". Then make a new type, maybe called ComponentFermion that takes an extra named argument ("Arg") which sets which name of QN to use. Our QN system will treat QN values with different names as separate, and separately conserve them. Make sure your Hamiltonian also conserves them too, otherwise you'll get an error.

Best,
Miles
commented by (360 points)
Hi again Miles,

Thanks for the answer.

About point 2): I have to say I do not know how to define this additional quantum numbers. Do we have any examples or info in the tutorials about it?

Rafael
commented by (70.1k points)
Hi Rafael,
Good question about a tutorial. We don't have one currently; it's on our list of things to add to the website.

In the meantime, almost better than a tutorial (in the sense of having a very closely related example to start from), your best bet would be to copy the file (from the library source code) itensor/itensor/mps/sites/fermion.h to a new file "component_fermion.h".

Then change the type name "FermionSite" everywhere in that file to "ComponentFermionSite" (if that's the new name you want to use). And similarly change the word "Fermion" to "ComponentFermion" in other cases. Also change the "#ifndef __ITENSOR_FERMION_H" lines and similar to use a different name of the "scope guard" macros there.

Now your simplest move will be to use the "SiteNumber" arg which gets passed to your ComponentFermionSite type to set the QN's used in the conserve_Nf == true case to either "Nfa", "Nfb", or "Nfc" depending on whether n%3 == 0, 1, or 2.

You can remove all other code (definitions of Index objects) that's not going to be used, such as the code related to only conserving parity, or not conserving any QNs at all.

Ok please take a try at that and we can discuss specific questions you may have as you go -

Miles
commented by (360 points)
Dear Miles,

Thank you so much for the detailed explanation! I will try to implement the model and see what I get. In case it works, I will probably share it here, I think this is something that may come up in future questions as well.

Best wishes,

Rafael
commented by (360 points)
Hello Miles,

I have implemented the code for a three-component multiband fermion model as we discussed here (I am still not enforcing the conservation in the number for each component, though).

The issue I have found is that, for some reason, running the dmrg calculations does not change the initial state at all. This is made evident by printing the fermion density in the system (see su3 code below).

This must be a pretty fundamental problem when implementing the hopping between every third site, because apparently there is no hopping at all. I tried playing around with bosonic operators and the fermionic strings in the Hamiltonian as well, but the result is the same.

On the other hand, if I revert the code to a single component case (getting rid of the "longer range" hopping) it seems to work fine (see spinless_ferimions code below).

Any ideas on what the issue might be?

-----------------------su3-----------------------

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

int main()
    {
    auto N = 27;//system size
    auto Nc = 3; //number of components
    float t = 1;
    float r = 1;
    float U = 0;
    
    println("number of sites ",N);
    println("hopping ",t);
  
    auto sites = Fermion(N);
    auto ampo = AutoMPO(sites);
    
    vector< pair <int,int> > u1prs;
    vector< pair <int,int> > t1prs;
    vector< pair <int,int> > t2prs;
    vector< pair <int,int> > t3prs;
    
    for(int i = 0; i < N; i=i+3)
        {
        u1prs.push_back(make_pair(i+0,i+1));
        u1prs.push_back(make_pair(i+0,i+2));
        u1prs.push_back(make_pair(i+1,i+2));
        }
        
    /*
    for(int i = 0; i < u1prs.size(); i++)
        {
        cout << u1prs[i].first << " ";     
        cout << u1prs[i].second << endl;
        }    
    */
    
    for(int i = 0; i < N-3; i=i+3)
        {
        t1prs.push_back(make_pair(i+0,i+3));
        t2prs.push_back(make_pair(i+1,i+4));
        t3prs.push_back(make_pair(i+2,i+5));
        }    
    /*
    for(int i = 0; i < t1prs.size(); i++)
        {
        cout << t1prs[i].first << " ";     
        cout << t1prs[i].second << endl;
        }   
    */
    
    for(int i = 0; i < u1prs.size(); i++)
        {
        int u1 = u1prs[i].first + 1;
        int u2 = u1prs[i].second + 1;
        println("Interactions pairs ", u1, " " ,u2);
        ampo += U,"N",u1,"N",u2;
        }
    
    for(int i = 0; i < t1prs.size(); i++)
        {
        int t1a = t1prs[i].first + 1;
        int t1b = t1prs[i].second + 1;
        println("Tunneling 1 ", t1a, " " ,t1b);
        ampo += -t,"Cdag",t1a,"C",t1b;
        ampo += -t,"Cdag",t1b,"C",t1a;
    }
   
    for(int i = 0; i < t2prs.size(); i++)
        {
        int t2a = t2prs[i].first + 1;
        int t2b = t2prs[i].second + 1;
        println("Tunneling 2 ", t2a, " " ,t2b);
        ampo += -t,"Cdag",t2a,"C",t2b;
        ampo += -t,"Cdag",t2b,"C",t2a;
        }
    
    for(int i = 0; i < t3prs.size(); i++)
        {
        int t3a = t3prs[i].first + 1;
        int t3b = t3prs[i].second + 1;
        println("Tunneling 3 ", t3a, " " ,t3b);
        ampo += -t,"Cdag",t3a,"C",t3b;
        ampo += -t,"Cdag",t3b,"C",t3a;
        }
    
    auto H = toMPO(ampo);
    
    auto state = InitState(sites);
    int p = 3 * Nc;
    for(int i = 1; i <= N; i++)
        {
        if(p > 0)
            {
            println("Singly occupying site ",i);
            state.set(i,"Occ");
            p -= 1;
            }
        else
            {
            println("Empty site ",i);
            state.set(i,"Emp");
            }
        }

    auto psi0 = MPS(state);  
    
    auto sweeps = Sweeps(10); //number of sweeps is 5
    sweeps.maxdim() = 10,20,100,100,200;
    sweeps.cutoff() = 1E-10;
    
    auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet"});
    println("Ground State Energy = ",energy);
    
    Vector d(N);
    
    for(int j = 1; j <= N; j++)
        {
        psi.position(j);
        d(j-1) = elt(dag(prime(psi(j),"Site"))*op(sites,"N",j)*psi(j));
        }
            
    println("Density:");
    for(int j = 0; j < d.size(); ++j)
        printfln("%d %.10f",1+j,d(j));
    println();

    return 0;
    }

-------------------------spinless_fermions-------------------------

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

int main()
    {
    auto N = 20;//system size
    auto t = 1;
    auto Nferm = 5;
    
    println("number of sites ",N);
    println("hopping ",t);
  
    auto sites = Fermion(N);
    auto ampo = AutoMPO(sites);
    
    vector< pair <int,int> > u1prs;
    vector< pair <int,int> > t1prs;
    
    for(int i = 0; i < N-1; i=i+1)
        {
        t1prs.push_back(make_pair(i,i+1));
        }    
        
    for(int i = 0; i < t1prs.size(); i++)
        {
        int t1a = t1prs[i].first + 1;
        int t1b = t1prs[i].second + 1;
        println("Tunneling pairs ", t1a, " " ,t1b);
        ampo += -t,"Cdag",t1a,"C",t1b;
        ampo += -t,"Cdag",t1b,"C",t1a;
    }
   
    auto H = toMPO(ampo);
    
    auto state = InitState(sites);
    int p = Nferm;
    for(int i = 1; i <= N; i=i+1)
        {
        if(p > 0)
            {
            println("Singly occupying site ",i);
            state.set(i,"Occ");
            p -= 1;
            }
        else
            {
            println("Empty site ",i);
            state.set(i,"Emp");
            }
        }

    auto psi0 = MPS(state);  
    
    auto sweeps = Sweeps(35);
    sweeps.maxdim() = 200,400,600,800;
    sweeps.cutoff() = 1E-10;
    
    auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet"});
    println("Ground State Energy = ",energy);
    
    Vector d(N);
    
    for(int j = 1; j <= N; j++)
        {
        psi.position(j);
        d(j-1) = elt(dag(prime(psi(j),"Site"))*op(sites,"N",j)*psi(j));
        }
        
    println("Density:");
    for(int j = 0; j < d.size(); ++j)
        printfln("%d %.10f",1+j,d(j));
    println();

    return 0;
    }
commented by (70.1k points)
Hi Rafael, I haven't had time to go through your code in detail yet, but a good way to imagine a DMRG calculation like this is to think of the particles as starting on whichever sites you put them on in the initial state, then "diffusing" around by following the hoppings you have defined. If certain sets of sites are disconnected from others, or only connected by a very weak hopping, then particles may never find their way to those sites. Other ways things can go wrong is if all of the particles start out piled up on one side of the system or only on one type of site (an "a" site versus a "b" or "c" site). Then it could take a very long time for the particles to make their way to other sites.

What I'm basically saying boils down to trying different initial states, as well as making sure there's nothing too pathological about your Hamiltonian (such as having disconnected sites).

Finally, something you haven't taken advantage of yet is the "noise term" feature of our DMRG code. If you add a line like:
sweeps.noise() = 1E-5,1E-6,1E-8,1E-10,1E-12;
then your code will run a bit slower but it will give you wavefunction a small, perturbative correction at each step that can greatly help convergence for complicated, QN-conserving cases like yours.

Please do these things then:
1. inspect your Hamiltonian to make sure it's ok
2. turn on the noise term
3. try some other initial states

One other good idea is doing lots of sweeps at smaller bond dimension before increasing the bond dimension. So you should do something like 10 sweeps at maxdim 20, then 10 at maxdim 40, then finally increase to maxdim's in the hundreds.

Best,
Miles
commented by (360 points)
Hi Miles,

Thanks again for the feedback.

I understand what you mean in terms of the "diffusion" of the initial state. It is instructive to notice that if only next-nearest neighbour hoppings are considered, there is no change to the initial state. That is true even for the "exthubbard" example contained in the sample folder.

One solution I found was to have a very small value for the nearest neighbour hopping which allows the particles to diffuse towards the ground state. But even better than that is turning on the noise term, as you said. Apparently just adding the values you suggested there is enough.

Now, what is still unclear to me is how to enforce the conservation of particle number for each component. My main issue with this is that I am adding an extra hopping term that only acts on components B and C, for example. This leads, in some cases, to a ground state where all particles localise in A sites, which is not what I want.

It still seems that the best approach would be to edit the change the "fermion.h" file, but it's not clear how to introduce in this file the distinction between particles and different sites. I assume that (besides changing the names of the type everywhere in the file) you were talking about editing the section below?

  if(conserve_Nf) //usual case
                {
                s = Index(QN({"Nf",0,-1}),1,
                          QN({"Nf",1,-1}),1,Out,ts);
                }
            else
                {
                s = Index(QN({"Pf",0,-2}),1,
                          QN({"Pf",1,-2}),1,Out,ts);
                }

Rafael
commented by (120 points)
Hi, I have closely related followup questions.  
 Is there any example in Julia for implementing multi-flavor Hubbard model?
commented by (70.1k points)
Hi Asad, we don't have any example of that that I know of. But here is how to do it: you just need to assign different sites to the flavors. So like if there are two flavors, then you could use odd numbered sites for flavor 1 and even for flavor 2. Then a first-neighbor hopping between sites of flavor 1 becomes a second-neighbor hopping under this mapping.

Does that give you what you are looking for?

Best regards,
Miles
commented by (120 points)
Thank you for your reply. I kind of did something like your solution. I label one flavor for site index 1 to N and another flavor for site index N+1 to 2N, and accordingly constructed different terms. I have one confusion. Say I have N sites, do I need different site indices objects for different flavors like
s1 = siteinds("S=1/2",N)
s2= siteinds("S=1/2",N)

or, is it fine if I just define one site indices object like this
s = siteinds("S=1/2",2*N)
commented by (70.1k points)
Hi Asad, would you be ok with asking this as a new question? Also it would be helpful if you posted which Hamiltonian you are wanting to simulate. We could discuss more on a new thread, for example it will be quite important not to organize the sites from 1..N then N+1,...2N but in a different pattern most likely.
commented by (120 points)
Hi Miles, thanks, I will post it as a new question.
Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.

Categories

...