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;
}