DMRG for obtaining multi excited states with conserved sz

Dear all,

I am currently working on applying DMRG to compute the ground-state energy and numerous excited state energies of a finite-size (N=6, 12) spin-1 chain subjected to a magnetic field (Bz), assuming conserved sz. Following this, I aim to determine the partition function to extract the free energy and magnetization at low temperatures. However, when utilizing the DMRG method, I consistently obtain zero magnetization for all values of the magnetic field.

Here’s the function that I’ve been using:

# Calculate free energy and magnetization
function calculate_FreeE_and_Mag(N, J, B_min, B_max, uB, sweeps, linkdim, cutoff, noise, T)
    sites = siteinds("S=1", N; conserve_sz=true)
    sB = Int((B_max - B_min) / uB)  # number of field steps
    h = zeros(sB + 1)                       # magnetic field
    F_values = zeros(sB + 1)          # free energy
    M_values = zeros(sB + 1)         # magnetization

    for (m_B, Bz) in enumerate(B_min:uB:B_max)
        h[m_B] = Bz

        statei = ["Emp" for n in 1:N]  # assuming conserved sz with permutation 3
        for i in 1:3:N
            statei[i] = "Up"
        end
        for i in 2:3:N
            statei[i] = "Z0"
        end
        for i in 3:3:N
            statei[i] = "Dn"
        end

        ψi = randomMPS(sites, statei, linkdim)
        H = H_XXXchain(N, J, J2z, Bz, sites)

        E0, ψ0 = dmrg(H, ψi, sweeps; write_when_maxdim_exceeds=1000)
        En = [E0]
        ψn = [ψ0]
        for i in 1:12  # number of exited states for N=12
            E, ψ = dmrg(H, ψn, ψi, sweeps; write_when_maxdim_exceeds=1000)
            @show i   E
            println()
            push!(En,E)
            push!(ψn,ψ)
        end
        # ------------------------------------------------------------------
        Z = sum(exp.(-En / T))   # partition function
        F_values[m_B] = -T * log(Z)   # free energy
        # Calculate magnetization using finite difference
        if m_B > 1
            M_values[m_B] = (1/N)*(F_values[m_B - 1] - F_values[m_B]) / uB  # magnetization
        end

        println("Bz = $Bz, Free Energy = $(F_values[m_B]), Magnetization = $(M_values[m_B])")
    end

    return F_values, M_values, h
end

However, when comparing this to the exact diagonalization (ED) method for N=12, the ED method yields accurate magnetization values for each magnetic field and temperature. I use following expression for ED method:

H_tensor = prod(H)
eigen_E, _ = eigen(H_tensor)
Z = sum(exp.(-eigen_E / T))
F_values[m_B] = -T * log(Z)
# Calculate magnetization using finite difference
if m_B > 1
    M_values[m_B] = (1/N)*(F_values[m_B - 1] - F_values[m_B]) / uB
end

Unfortunately, in my case, ED does not work for N>=12. In fact, I get following error:

Contraction resulted in ITensor with 15 indices, which is greater
            than or equal to the ITensor order warning threshold 14.
            You can modify the threshold with macros like `@set_warn_order N`,
            `@reset_warn_order`, and `@disable_warn_order` or functions like
            `ITensors.set_warn_order(N::Int)`, `ITensors.reset_warn_order()`, and
            `ITensors.disable_warn_order()`.
Stacktrace:
  [1] _contract(A::ITensor, B::ITensor)
    @ ITensors ~/.julia/packages/ITensors/Gf9aD/src/tensor_operations/tensor_algebra.jl:20
  [2] contract(A::ITensor, B::ITensor)
    @ ITensors ~/.julia/packages/ITensors/Gf9aD/src/tensor_operations/tensor_algebra.jl:104
  [3] *
    @ ~/.julia/packages/ITensors/Gf9aD/src/tensor_operations/tensor_algebra.jl:91 [inlined]
  [4] mul_prod
    @ ./reduce.jl:35 [inlined]
  [5] BottomRF
    @ ./reduce.jl:86 [inlined]
  [6] _foldl_impl
    @ ./reduce.jl:62 [inlined]
  [7] foldl_impl(op::Base.BottomRF{typeof(Base.add_sum)}, nt::Base._InitialValue, itr::ITensor)
    @ Base ./reduce.jl:48 [inlined]
  [8] mapfoldl_impl
    @ ./reduce.jl:44 [inlined]
  [9] mapfoldl
    @ ./reduce.jl:175 [inlined]
 [10] mapreduce
    @ ./reduce.jl:307 [inlined]
 [11] prod
    @ ./reduce.jl:620 [inlined]

Despite having 64 gigabytes of RAM, assuming conserved sz, and implementing the term ITensors.disable_warn_order(), ED for N=12 spin-1 particles unexpectedly halts due to memory limitations.

My guess is that the DMRG should also provide accurate results when sz is conserved. It’s worth noting that when I consider non-conserved sz, the DMRG works and provides accurate magnetization for any magnetic field and temperature. However, when conserved sz is considered, the computation time significantly decreases.

I’ve also tried with the predefined code purification.jl, but unfortunately, I couldn’t obtain accurate magnetization values for non-zero magnetic fields. Despite this, using the thermal partition function, I can still obtain other thermodynamic quantities.

If you have any insights or suggestions on how to address the issues for ED and DMRG and improve the accuracy of the magnetization results in the DMRG method, please let me know.

Thanks