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 BoseHubbard chain in the Mottinsulating 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 unitfilling (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:(L1)
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, 1E8)
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_bar1.)*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