I am trying to reproduce the result for ferromagnetism in 1d chain using DMRG particularly from this paper DMRG study of ferromagnetism in a one-dimensional Hubbard model | SpringerLink. I am not able to figure out the problem with my code. It is giving me spin_up number same as spin_down instead of giving me a ferromagnetic ground state
‘’’
using ITensors
using LinearAlgebra
using Plots
a=Float64
b=Float64
let
L=20
U=5
h=0
N = 40
sites = siteinds(“Electron”, N;conserve_nf=true)
@show sites[1]
t =-1
tt=0
mu=0
ampo = OpSum()
for j=1:N-1
ampo += t,“Cdagup”,j,“Cup”,j+1
ampo += t,“Cdagup”,j+1,“Cup”,j
ampo += t,“Cdagdn”,j,“Cdn”,j+1
ampo += t,“Cdagdn”,j+1,“Cdn”,j
ampo += U,“Nup”,j,“Ndn”,j
ampo += mu,“Ntot”,j
ampo += h,“Sz”,j
end
for j=1:N-2
ampo += tt,“Cdagup”,j,“Cup”,j+2
ampo += tt,“Cdagup”,j+2,“Cup”,j
ampo += tt,“Cdagdn”,j,“Cdn”,j+2
ampo += tt,“Cdagdn”,j+2,“Cdn”,j
end
ampo += h,“Sz”,N
ampo += mu,“Ntot”,N
ampo += U,“Nup”,N,“Ndn”,N
state1 = [isodd(n) ? “Emp” : “Up” for n=1:N]
H = MPO(ampo, sites);
sweeps = Sweeps(100);
maxdim!(sweeps, 64,64,100,100,128,128,160,160,200,200,256,256,300,300,400,400,400);
setcutoff!(sweeps,1E-15);
setnoise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)
psiii = randomMPS(sites,state1,8);
energy, psi = dmrg(H, psiii, sweeps);
expNup=expect(psi,“Nup”);
expNdn=expect(psi,“Ndn”);
nt=expect(psi,“Ntot”)
println(“$nt”)
Szexp=expect(psi,“Sz”)
nbar=sum(expNup)/N;
ntot=sum(expNdn)/N;
println(“$nbar”)
println(“$ntot”)
println(“$Szexp”)
zz=correlation_matrix(psi,“Sz”,“Sz”);
z1=diag(zz)
sss=size(zz)
println(“$z1”)
‘’’
could you please help me figure out the mistake?