Hello iTensor community,
I just started to use iTensor few days ago and I have several questions about the DMRG part (Julia version). My goal is to recover easy results concerning the Bose-Hubbard chain in the Mott-insulating phase. More precisely, my goal is to find after running the DMRG algorithm the correct ground state corresponding to a Mott state starting from an initial guess (a product state) permitting to fix the total number of bosons on the lattice. However, I have several problems (see Julia code at the end of the message):
-
Firstly, my code does not conserve the total number of bosonic particles on the lattice after running the DMRG function. I do not understand since I defined initially a product state permitting to fix the number of particles and besides the Hamiltonian conserves the total number of particles.
To solve the problem I tried to add:conserve_number = true
or ’conserve_qns = true
leading for example to:sites = siteinds("Boson", L; conserve_number = true, dim = n_bar*L+1)
but I obtained an error during the compilation in both cases. Do you have an idea concerning the origin of the mistake? -
Secondly and maybe if the first point explained above is solved will solve the following questions: I do not understand why I should consider such a high value of
dim
corresponding to the maximum number of bosons on a lattice site (if I am correct). To get the correct ground state, It should be enough to considerdim = n_bar+1
since the ratio U/J is much larger (U/J = 20 here) than the critical value at unit-filling (n_bar = 1) given by (U/J)c = 3.3. Same remark for the number of sweeps in order to reach the convergence. -
I do not know and I did not find in the documentation the meaning of
noise
. Besides, I do not understand why the number of values ofmaxdim
andnoise
do not match the number of sweeps.
I hope that these questions are not trivial and will be helpful for the iTensor community.
Best regards,
Julien DESPRES
using ITensors
let
L = 20
n_bar = 1
J = 1.
U = 20.
sites = siteinds("Boson", L; dim = n_bar*L+1)
os = OpSum()
for j = 1:(L-1)
os += - (U/2.), "N", j
os += (U/2.), "N", j, "N", j
os += -J, "Adag", j, "A", j+1
os += -J, "A", j, "Adag", j+1
end
os += - (U/2.), "N", L
os += (U/2), "N", L, "N", L
# periodic boundary conditions
# os += -J, "A", 1, "Adag", L
# os += -J, "Adag", 1, "A", L
# println(os)
H = MPO(os, sites)
states = ["2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0"]
psi0 = MPS(sites, states)
# println(psi0)
dens = expect(psi0, "N")
N_tot_init = 0.
for (j, dn) in enumerate(dens)
println("$j $dn")
N_tot_init += dn
end
println("N_tot_init = $N_tot_init")
sweeps = Sweeps(200)
setmaxdim!(sweeps, 10, 20, 40, 80, 160)
setcutoff!(sweeps, 1E-8)
setnoise!(sweeps, 0.0, 0.0, 0.0, 0.0)
energy, psi = dmrg(H, psi0, sweeps)
println("E_GS = $energy")
energy_th = (U/2.)*n_bar*(n_bar-1.)*L
println("E_GS_th_U_inf = $energy_th")
dens = expect(psi, "N")
N_tot_final = 0.
for (j, dn) in enumerate(dens)
println("$j $dn")
N_tot_final += dn
end
println("N_tot_final = $N_tot_final")
return
end