How to construct phase diagram of bose-hubbard model using Julia?

ITensor Julia Questions
1

Hello everyone!

I’d like to ask a question about how to set the total number of particles in bose-hubbard model “without controlling the particle number of every site”.

Its Hamiltonian is
\hat{H}=-t\sum*{i=1}^L(\hat{b}*i^\dagger \hat{b}*{i+1}+H.c.)+\frac{U}{2}\sum*{i=1}^L\hat{n}*i(\hat{n}*i-1)-\mu\sum*{i=1}^L\hat{n}*i

I want to calculate its phase diagram through its ground-state energy difference when one particle was added or removed.
\mu(N)_+=E(N+1)-E(N)
\mu(N)_-=E(N)-E(N-1)

Here is my code by Julia version. I want to know how my code runs correctly and control the particle number on every site but I just want to control the total number of particles, not every site.

using ITensors

J = parse(Float64, ARGS[1])
miu = parse(Float64, ARGS[2])
U = parse(Float64, ARGS[3])
V = parse(Float64, ARGS[4])
L = parse(Int, ARGS[5])

sites = Boson(L, "MaxOcc=", 5, "ConserveQNs", true, "ConserveNb", false)

sweeps = Sweeps(8)
sweeps.maxdim() .= [40, 80, 400, 400, 800, 800, 1000, 1000]
sweeps.cutoff() = 1E-16

ampo = AutoMPO(sites)
for i in 1:L-1
    add!(ampo, -miu - U/2, "N", i)
    add!(ampo, U/2, "N", i, "N", i)
    add!(ampo, -J, "A", i, "Adag", i + 1)
    add!(ampo, -J, "Adag", i, "A", i + 1)
    add!(ampo, V, "N", i, "N", i + 1)
end
add!(ampo, -miu - U/2, "N", L)
add!(ampo, U/2, "N", L, "N", L)

H = toMPO(ampo)

state = InitState(sites)
for i in 1:L
    setindex!(state, "3", i)
end

PrintData(sweeps)

psi0 = randomMPS(state)
energy, psi = dmrg(H, psi0, sweeps, "Quiet=", true)

file_name1 = "your_file1.txt"
open(file_name1, "w") do fp
    println(fp, "$miu,$energy")
end

file_name2 = "your_file2.txt"
open(file_name2, "w") do fp
    println(fp, "i,j,lgj,corf")
    
    Numall = 0
    for i in 1:L
        position(psi, i)

        ket = psi[i]
        bra = dag(prime(ket, "Site"))

        Numop = op(sites, "N", i)
        Numi = elt(bra * Numop * ket)
        Numall += Numi
        println(fp, "$i,%.12f", Numi)
    end
    
    println(fp, "i,j,lgj,corf")
    
    i = 1
    for j in 2:L
        position(psi, i)

        corfop_i = op(sites, "Adag", i)
        corfop_j = op(sites, "A", j)

        psidag = dag(psi)
        psidag.prime("Link")
        l_i = leftLinkIndex(psi, i)
        C = prime(psi[i], l_i) * corfop_i
        C *= prime(psidag[i], "Site")

        for k in i+1:j-1
            C *= psi[k]
            C *= psidag[k]
        end

        lj = rightLinkIndex(psi, j)
        C *= prime(psi[j], lj) * corfop_j
        C *= prime(psidag[j], "Site")
        corf_1j = elt(C)
        println(fp, "$i,$j,$(log10(j)),$corf_1j")
    end
end

return 0

Looking forward to your help.

2

If you use quantum number conservation in ITensor, the behavior will actually be to conserve the total number of particles. It will not conserve a certain number of particles “on every site”, as separate conservation laws.

In the Julia version, you can construct a set of site indices with Boson properties as follows:

s = siteinds("Boson",N; conserve_qns=true, dim=4)

Above I have set the site dimension to 4, but you can set it to other values to use a larger (truncated) bosonic Hilbert space. Setting conserve_qns=true will add QN labels to each site, but then the thing that is conserved is the total of all these labels, not the individual labels.

The way that you specify the number of particles you want there to be in the system is through the initial state you construct as input to dmrg. So when you create the vector state the number of particles will be set as the total across the whole initial state.

Finally, to get the energies that you want such as E(N+1), you can just do separate DMRG calculations, one for each total number of particles. Then you can subtract the various energies to get the chemical potentials such as \mu_+ that you are aiming to compute.

Please see the fully working example of how to set up a quantum number conserving DMRG calculation in Julia here:
https://itensor.github.io/ITensors.jl/dev/tutorials/QN_DMRG.html

3

I’m sorry my code still doesn’t work, If I wanna picture the phase diagram using t_miu_superfluid order parameters as x, y, and z axes, What should I define the order parameter using ITensors. but for traditional DMRG, we can only use /mu_{+} and /mu_{-} to picture the phase diagram of the Mott to SF. If this, how do separate DMRG calculations, one for each total number of particles?

4

I think you are really asking some separate questions here:

  1. When you say your code doesn’t work, are you getting an error message? If so, what is it? Or do you mean the output is not what you expected or hoped to get? Please give enough details so that we can help you.

  2. About order parameters, there are various ones that people may try to use to characterize phases of matter. Most of them can be measured using tensor networks. What is the order parameter you are hoping to use? It is not a question one can answer from the point of view of what ITensor or DMRG can do, but is primarily a question about physics or your own personal research goals.

  3. To do separate DMRG calculations for different numbers of particles, call the dmrg routine each time, providing a different initial state which has the number of particles that you want. So for example if you want to compute energies for the cases of (N+1), N, and (N-1) particles, you would do three separate DMRG calculations. Please read the ITensor example/tutorial page I linked for more background on how our QN-conserving DMRG system works. (Basically whatever total initial QN’s your initial state has are always preserved throughout DMRG. In your case it would be particle number. You can verify this by computing the total density using code like dens = expect(psi,"N"); @show sum(dens); before and after DMRG runs to verify what the particle number is and that it remains the same before and after.

5

The error messages are as follows:
julia “e:\exercise\Julia_learning\Bose-Hubbard.jl”
ERROR: LoadError: BoundsError: attempt to access 0-element Vector{String} at index [1]
Stacktrace:
[1] getindex(A::Vector{String}, i1::Int64)
@ Base .\essentials.jl:13
[2] top-level scope
@ e:\exercise\Julia_learning\Bose-Hubbard.jl:3
in expression starting at e:\exercise\Julia_learning\Bose-Hubbard.jl:3
I am beginning to learn ITensors and DMRG, It’s my honor to get your reply. Next, I will take your advice to learn. I sincerely hope that I can get your answers when I encounter questions in the future.

6

If you read this kind of error message, it already has a lot of information. For example, it is telling you that the error came from line 3 of your program. Then, the error is that a vector (or array) was accessed at an entry that was larger (“out of bounds”) than its size. So please look at or near line 3 of your code and see what vector is being used there, and I would recommend printing out the size or length of this vector to see why the out of bounds access happened.