Dear there,
I’m playing with the simple 2-leg Hubbard model. I attach my short code below. And I got a problem that is when I run the half filling case (I have 16x2 sites, so 32 particles in total), there’s no problem, but when I trying case like 1 electron per two sites (16 particles total), all my sweeps show maxlinkdim=1, even though I’ve set both max and min dim to be 256:
setmaxdim!(sweeps_tem,256) , setmindim!(sweeps_tem,256)
I’ve also attached this error results. Could you help to see what I’m doing wrong with this standard Hubbard model?
Thank you so much!!
flux(psi0) = QN(("Nf",16,-1),("Sz",0))
sweeps_tem = Sweeps
1 cutoff=1.0E-16, maxdim=512, mindim=512, noise=0.0E+00
2 cutoff=1.0E-16, maxdim=512, mindim=512, noise=0.0E+00
3 cutoff=1.0E-16, maxdim=512, mindim=512, noise=0.0E+00
4 cutoff=1.0E-16, maxdim=512, mindim=512, noise=0.0E+00
5 cutoff=1.0E-16, maxdim=512, mindim=512, noise=0.0E+00
6 cutoff=1.0E-16, maxdim=512, mindim=512, noise=0.0E+00
After sweep 1 energy=0.0 maxlinkdim=1 maxerr=0.00E+00 time=19.378
After sweep 2 energy=0.0 maxlinkdim=1 maxerr=0.00E+00 time=0.017
After sweep 3 energy=0.0 maxlinkdim=1 maxerr=0.00E+00 time=0.026
After sweep 4 energy=0.0 maxlinkdim=1 maxerr=0.00E+00 time=0.017
After sweep 5 energy=0.0 maxlinkdim=1 maxerr=0.00E+00 time=0.025
After sweep 6 energy=0.0 maxlinkdim=1 maxerr=0.00E+00 time=0.016
using ITensors
using HDF5
let
N = 16
Npart = 16
t = 1.0
U = 1.0
sites = Vector{Index}(undef,2N)
for j=1:2N
sites[j] = siteind("Electron";addtags="n=$j",conserve_qns=true)
end
ampo = OpSum()
# x direction:
for n in 0:1
for i in n+1:2:(n+1+2*(N-2))
ampo += -t, "Cdagup", i, "Cup", i + 2
ampo += -t, "Cdagup", i + 2, "Cup", i
ampo += -t, "Cdagdn", i, "Cdn", i + 2
ampo += -t, "Cdagdn", i + 2, "Cdn", i
end
end
# y direction:
for i in 1:2:31
ampo += -t, "Cdagup", i, "Cup", i + 1
ampo += -t, "Cdagup", i + 1, "Cup", i
ampo += -t, "Cdagdn", i, "Cdn", i + 1
ampo += -t, "Cdagdn", i + 1, "Cdn", i
end
# Hubbard U interaction
for i in 1:2N
ampo += U, "Nup", i, "Ndn", i
end
H = MPO(ampo, sites)
state = ["Emp" for n in 1:2N]
p = Npart
count=1
# Half filling case (no problem)
for i in 1:2N
state[i]= (isodd(count) ? "Up" : "Dn")
count += 1
end
# 1 electrons per two sites (get stuck at maxlinkdim=1)
for i in 1:4:31
println("$i")
state[i]= (isodd(count) ? "Up" : "Dn")
count += 1
end
for i in 2:4:32
println("$i")
state[i]= (isodd(count) ? "Up" : "Dn")
count += 1
end
psi0 = productMPS(sites, state)
@show flux(psi0)
# Start DMRG calculation:
sweeps_tem= Sweeps(10)
bond=256
setmaxdim!(sweeps_tem,bond)
setmindim!(sweeps_tem,bond)
@show sweeps_tem
energy, psi_tem = dmrg(H, psi0, sweeps_tem; eigsolve_tol=1e-6)
end