DMRG - Conservation of the total number of particles for the Bose-Hubbard chain

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 consider dim = 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 of maxdim and noise 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
sweeps = Sweeps(10)
setmaxdim!(sweeps,1,2,3)
julia> @show sweeps
sweeps = Sweeps
1 cutoff=1.0E-16, maxdim=1, mindim=1, noise=0.0E+00
2 cutoff=1.0E-16, maxdim=2, mindim=1, noise=0.0E+00
3 cutoff=1.0E-16, maxdim=3, mindim=1, noise=0.0E+00
4 cutoff=1.0E-16, maxdim=3, mindim=1, noise=0.0E+00
5 cutoff=1.0E-16, maxdim=3, mindim=1, noise=0.0E+00
6 cutoff=1.0E-16, maxdim=3, mindim=1, noise=0.0E+00
7 cutoff=1.0E-16, maxdim=3, mindim=1, noise=0.0E+00
8 cutoff=1.0E-16, maxdim=3, mindim=1, noise=0.0E+00
9 cutoff=1.0E-16, maxdim=3, mindim=1, noise=0.0E+00
10 cutoff=1.0E-16, maxdim=3, mindim=1, noise=0.0E+00
1 Like

Also note that you no longer are required to make a Sweeps object to pass into dmrg and can now just pass numbers and arrays like

nsweeps = 5
maxdim = [10,20,100,100,200]
cutoff = [1E-10]

where as Ryan said, if the arrays are shorter than the number of sweeps, the last value in the array is implicitly repeated (reused) for the remaining sweeps.

Then you can pass these numbers and arrays as keyword arguments into dmrg like this:

energy, psi = dmrg(H,psi0; nsweeps, maxdim, cutoff)

Please see these examples which use this newer style and also the second example shows how to upgrade the first one to use quantum number conservation:

Hello Ryan, Miles,

Upgrading my Ubuntu version and reinstalling ITensors has solved the problem. Thank you very much for your reply and the information.

Best,
Julien