Hubbard model in magnetic field

+1 vote

Hello,

I was trying to simulate 1D Hubbard model ground state in presence of a magnetic field. When I vary the strength of the magnetic field, the ground state energy was changing, but the 1-site reduced density matrix does not change, which should not be the case.

I am writing the code below. Any suggestion will be appreciated.

auto t1 = 1;
auto U1 = 0;

auto sites = Electron(N);
auto ampo = AutoMPO(sites);

for(double mu=-4; mu<=4; mu=mu+0.5)
{

for(int i = 1; i <= N; ++i)
{
ampo += U1,"Nupdn",i;
}
for(int i = 1; i <= N; ++i)
{
ampo += -mu,"Sz",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 psi0 = MPS(state);
Print(totalQN(psi0);
auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet",quiet});

mat rho_1 = RDM_site(psi, N/2, 4);
cout << rho_1 << endl;

}


where RDM_site is the function for calculating the 1-site reduced density matrix.

Hi Streetmanda,
By default, the indices created in an Electron site set keep track of the Sz quantum number. So then an MPS made from these indices will be a state of well-defined Sz that will not be changed, since your Hamiltonian conserves total Sz. Throughout the calculation, the total Sz will remain at the value of the input state psi0 used to initialize the calculation. This is how DMRG calculations with conserved total quantum numbers are carried out in ITensor.

So if you'd like the magnetic field to be able to alter the magnetization, you'll need to turn off Sz conservation, like this:

auto sites = Electron(N,{"ConserveSz=",false});


when making your site set sites.

For a list of the available options to the Electron site set you can see the documentation here:
http://itensor.org/docs.cgi?vers=cppv3&page=classes/electron

Best,
Miles

commented by (190 points)
Thank you very much Miles for your answer. I actually tried that before posting the question. But even when using this, the reduced density matrix does not change for certain segments. For example, when the magnetic filed strength "mu" is changed from -1 to +0.8 in small  steps, the reduced density matrix remains same in all intermediate steps, which should not be the case. The energy is well converged in those steps. I will check my code again, but any insight will be helpful.
commented by (46.8k points)
Thanks for pointing this out. However, it's hard to answer this part of your question with the information provided. Are you sure that the single-site reduced density matrix does not change? (I believe you that is doesn't, but it's kind of puzzling that you could be getting the correct energy yet the properties of the wavefunction aren't changing.) Is your code to compute this RDM correct? Does that code call psi.position(j) to change the orthogonality center of the MPS? And you are saying it doesn't change even with total-Sz conservation turned off? Can you please test on a very small system (4 sites, say) against the exact results to test each part of your code?

Best,
Miles
commented by (190 points)
I am attaching the part of my code that I have changed with {"ConserveSz=",false}. I am also attaching the function I used for calculating the 1-site reduced density matrix. Yes I am getting a different energy for different field strengths, but the 1-RDM as well as the difference between up-spin and down-spin in each site remains unchanged for particular segments.

mat RDM_site(MPS psi, int site1, int d) {

psi.position(site1);  // 'gauge' the MPS to site site1
auto psidag = dag(psi);
auto rho = psi(site1)*prime(psidag(site1),"Site");
// PrintData(rho);

mat rho_1(d,d);
for(int i = 1; i <= d; ++i){
for(int j = 1; j <= d; ++j){
rho_1(i-1,j-1)=rho.real(i,j);
}
}

return rho_1;
}

The site command is now,

auto sites = Electron(N,{"ConserveSz=",false});
auto ampo = AutoMPO(sites);

and I am using the following steps to check the difference between up-spin and down-spin at each site:

Vector upd(N),dnd(N);
for(int j = 1; j <= N; ++j)
{
psi.position(j);
upd(j-1) = elt(dag(prime(psi(j),"Site"))*op(sites,"Nup",j)*psi(j));
dnd(j-1) = elt(dag(prime(psi(j),"Site"))*op(sites,"Ndn",j)*psi(j));
}

println("Density Difference:");
for(int j = 0; j < N; ++j)
printfln("%d %.10f",1+j,(upd(j)-dnd(j)));
println();

I will try to match it with exact results for small number of sites. I tried changing the initial state but that does not make any difference.
commented by (46.8k points)
Hi Prakash,
It's still a bit hard to say, but my main suspicion is how you are doing multiple DMRG calculations within a loop. I have done this myself in the past, and have seen other users on here do it, and it almost never goes well. The reason is that DMRG calculations take some care to converge for tougher systems, such as the Hubbard model, as you probably know. For the Hubbard model especially, DMRG can be very sensitive to the choice of initial state, and it's often quite important to both turn on the "noise" parameter (noise term) during certain sweeps as well as doing many sweeps. Thus it's best to study one particular Hamiltonian (value of mu) at a time to ensure convergence, and go through values of mu one at a time after that.

In short, I believe a likely reason for what you're seeing is that your DMRG calculations are not well converged here, or are getting stuck and not yet reaching the ground state.

Best,
Miles