mean of normalized engeries gaps

[Veneme]

I am a beginner and I am using ITensors.jl to calculate the many-body localization properties of the gradient field of 1D Bose-Hubbard model. However, I have encountered a problem: the results of my DMRG calculation do not match with the results obtained using Quspin benchmarks. But I seem to feel that there is no problem with my program. Here is the code I am using with ITensors as well as the code and the figure from Quspin.

using ITensors
using Statistics

L = 10
Np = 5
# Create five bosonic sites
sites = siteinds("Boson", L; conserve_qns=true, dim=4)

# Define other Hamiltonian parameters
t = 1
U = 7 # interaction
mu = 0 # chemical potential

# Define the range of g values to consider
gs = collect(0.1:0.1:2.1)

# Loop over g values and compute ground state energies
energies = Float64[]
diff_energies = Float64[]
r_spacing = Float64[]
for g in gs
    tg =0.5 * g * (1 / sinh(g))*(-t)
    # Compute the modified hopping term for this g value
    os = OpSum()
    for j in 1:(L - 1)
        os += g*j -mu - U*0.5, "N", j
        os += U*0.5, "N", j, "N", j
        os += tg, "Adag", j+1, "A", j
        os += tg, "Adag", j, "A", j + 1
    end
    os += g*L-mu - U*0.5, "N", L #OBC
    os += U*0.5, "N", L, "N", L
    # os += tg, "A", 1, "Adag", L #PBC
    # os += tg, "Adag", 1, "A", L
    H = MPO(os, sites)

    # Create a random initial state
    states = ["0" for j=1:L]
    for j in 1:Np
        states[j] = "1"
    end
    psi0 = randomMPS(sites, states)

    # Compute the ground state energy
    sweeps = Sweeps(15)
    setmaxdim!(sweeps, 10,20,40,60,100,140,200,250,300,350,400,450,500,550,600)
    setcutoff!(sweeps, 1E-13)
    setnoise!(sweeps, 1E-7, 1E-10, 1E-5, 0, 1E-8, 1E-7)
    
    energy, psi = dmrg(H, psi0, sweeps)
    
    # Append the energy to the list of energies
    push!(energies, energy)
end

###########
# Define the diff_energies array
diff_energies = [energies[i+1] - energies[i] for i in 1:length(energies)-1]

# Define the normal_r array
normal_r = [min(diff_energies[i-1], diff_energies[i]) / max(diff_energies[i-1], diff_energies[i]) for i in 2:length(diff_energies)]

# Calculate the mean of normal_r
r_mean = mean(normal_r)

println("$(r_mean)")
# Write g and r_mean arrays to a.txt file
open("b.txt", "w") do f
    for (g, r) in zip(gs[1:end-1], normal_r)
        println(f, "$g $r")
    end
end

Here is the DMRG calculation result. The second line represents the mean of normalized energy gaps, which seems to have no significant changes.

Quspin code：

from quspin.operators import hamiltonian
from quspin.basis import boson_basis_1d
import numpy as np
import matplotlib.pyplot as plt

# define model parameters
t = 1.0 # hopping
alpha = 0.0
beta = 0.0
gamma_vals = np.arange(0.1, 5.0, 0.1) # range of gamma values
U = 0.0 # interaction
mu = 2.71 # chemical potential

# initialize arrays to store properties
mean_normalized_gaps = []
N_vals = range(10, 15)

# loop over N values
for N in N_vals:
    # initialize arrays to store properties
    energies = []
    normalized_gaps = []
    
    # loop over gamma values
    for gamma in gamma_vals:
        # define site-coupling lists
        def g(i):
            return alpha*i**2 + beta*i + gamma
        tg = 0.5*g(0)*np.reciprocal(np.sinh(g(0)))*(-t) # new hopping parameter
        hop=[[tg,i,(i+1)] for i in range(N-1)] #OBC
        for i in range(N-1):
            tg = 0.5*g(i+1)*np.reciprocal(np.sinh(g(i)))*(-t)
            hop[i][0] = tg
        interact=[[0.5*U,i,i] for i in range(N)] # U/2 \sum_j n_j n_j
        pot=[[g(i)*i-mu-0.5*U,i] for i in range(N)] # -(\mu + U/2) \sum_j j_n
        # define static and dynamic lists
        static=[['+-',hop],['-+',hop],['n',pot],['nn',interact]]
        dynamic=[]
        # build Hamiltonian
        basis = boson_basis_1d(N,Nb=N) # full boson basis
        H = hamiltonian(static,dynamic,basis=basis,dtype=np.float64)
        # calculate eigensystem
        E, V = H.eigh()
        energies.append(E) # store eigenenergies
        gaps = np.diff(E)
        # calculate normalized energy gaps
        normalized_gap = np.min([gaps, np.roll(gaps,-1)], axis=0) / np.max([gaps, np.roll(gaps,-1)], axis=0)
        normalized_gaps.append(normalized_gap[:-1]) # exclude last element
    
    # calculate mean of normalized energy gaps
    mean_normalized_gaps.append(np.mean(normalized_gaps, axis=1))
    # Save mean_normalized_gaps and gamma_vals to file
    np.savetxt("level-spacing.txt", np.column_stack((gamma_vals, np.array(mean_normalized_gaps).T)), header="gamma   normalized energy gaps", fmt='%-10.5f')
    # plot results
    markers = ['o', 's', 'D', '*', '^', 'v']
    colors = ['r', 'g', 'b', 'm', 'c', 'k']
    for i in range(len(N_vals)):
    plt.plot(gamma_vals, mean_normalized_gaps[i], marker=markers[i], color=colors[i], label='N='+str(N_vals[i]))
    plt.xlabel(r'$\gamma$')
    plt.ylabel('means_of_$r$')
    plt.legend()
    plt.show()

Thank you very much for taking the time to read this！

[miles]

Glad you’re trying out ITensor. Regarding questions of best practices for using numerical methods or how to ensure that ITensor DMRG calculations are converged, I’d recommend the following FAQ section of our docs:
https://itensor.github.io/ITensors.jl/stable/faq/DMRG.html#Ensuring-a-DMRG-calculation-is-converged

In general I would strongly discourage putting a series of DMRG calculations inside of a loop over parameters. It’s important to ensure that DMRG is well converged for each Hamiltonian separately, by varying the number of sweeps and other parameters such as the maxdim and cutoff. I’d encourage you also, for a given Hamiltonian, to plot the properties of the state in real space to see if they look reasonable (smooth, left-right symmetric, etc.).

Also in general with DMRG, the first thing to try is just to do more sweeps. (Similar to how in Monte Carlo one would collect more samples). Hope that helps –

[Veneme]

Hello Miles! Thank you very much for your answers and guidance. They have been extremely helpful to me, and I am currently preparing to learn how to solve these problems.