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.