Conserving total quantum number SzT in spinful fermionic ladder

ITensor Julia Questions
1

Hi everyone,

My goal is to investigate the ground-state energy of a spinful Hubbard ladder versus the total SzT per site (E_0 vs. SzT).

In the case of spin-1/2 Heisenberg models, I only need to consider the total magnetization per site as a conserved quantity that ranges from zero to N/2. To achieve this, I define the lattice sites as follows:

sites = siteinds("S=1/2",N; conserve_sz=true)

I initialize the state statei as:

statei = ["Dn" for n in 1:N]       # all spins down
for j in 1:Int(N/2)                     # producing antiferromagnetic state  
     statei[j] = "Up"

This setup creates a fully antiferromagnetic state with SzT = 0. To calculate the quantity Szi, I define the following function:

# Sz(i) operator
function Szi_op(i,sites)
    ampo = AutoMPO()
    ampo += 1.,"Sz",i
    Szi = MPO(ampo,sites)
    return Szi
end

The main part of the dmrg code for the Heisenberg spin model can be summarized as:

spin_up = 0
for i in Int(N/2)+1: N
    spin_up += 1
    for j in i                # in each step a Dn in statei changes to Up
        statei[j] = "Up"
    end
    ψi = randomMPS(sites,statei,linkdim)
    Szi = [Szi_op(szi,sites) for szi in 1:N]
    H = H_Spin_Model(N,J,sites)     # Typical Hamiltonian's Function 
    E0,ψ0 = dmrg(H,ψi,sweeps)
    println()
    E_GS[spin_up] = E0
    # on-site total magnetization of the ground state ψ0 which is a conserved quantity
    for k in 1:N
        SzT[m_up] += inner(ψ0,Szi[k],ψ0)
    end
end

In this setup, the dmrg code runs successfully for Heisenberg spin models, producing the ground-state energy and on-site total magnetization for a sequence of states transitioning from fully antiferromagnetic to fully polarized (first for loop).

Now, I am attempting to extend the same methodology to the spinful Hubbard model. However, when I apply the previously described code, I am encountering an issue where the calculated ground-state energy turns out to be always zero.
In fact, I observe an array containing N/2 zeros representing the ground-state energy.

I am seeking guidance on how to address this issue and obtain accurate results for ground-state energy and the on-site magnetization of the spinful Hubbard model. Any insights or suggestions would be greatly appreciated.

2

It sounds like mainly you are having issues getting DMRG to converge. I would recommend the following steps to start with:

  1. Please try all of the recommendations in the FAQ linked below for ensuring DMRG convergence:
    https://itensor.github.io/ITensors.jl/dev/faq/DMRG.html

  2. I would strongly recommend against running multiple DMRG calculations inside of a for loop over different systems (different Hamiltonians and/or initial states). Instead, please study each system or initial state carefully and individually, by doing things like varying the number of sweeps or making plots of the results (e.g. plots of the electron density) to convince yourself that everything is converged, before going to the next system.

  3. It could also be good in your case to compute the exact solution on small systems within various spin or particle number sectors to check what the exact energy is and then compare this to your DMRG.

3

Dear miles,

When I use the exact solution for spin-1 chain with N=6, the ED method yields accurate magnetization values for each magnetic field and temperature. I use following code:

# Calculate free energy and magnetization
function calculate_free_energy_and_magnetization(N, J1, J2xy, J2z, B_min, B_max, uB, linkdim, T)
    sites = siteinds("S=1", N; conserve_qns=true, conserve_sz=true)
    # ITensors.disable_warn_order()
    sB = Int((B_max - B_min) / uB)
    h = zeros(sB + 1)
    F_values = zeros(sB + 1)
    M_values = zeros(sB + 1)

    for (m_B, Bz) in enumerate(B_min:uB:B_max)
       h[m_B] = Bz
        ψi = randomMPS(sites, linkdim)
        H = H_DiamondS1(N, J1, J2xy, J2z, Bz, sites)
        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

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

    return F_values, M_values, h
end

But when I increase the number of spins up to N=12, 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]
 [12] calculate_free_energy_and_magnetization(N::Int64, J1::Int64, J2xy::Float64, J2z::Float64, B_min::Int64, B_max::Int64, uB::Float64, linkdim::Int64, T::Float64)
    @ Main ~/ITensor Projects/FiniteT/ED_localSz_Sz_B_T.jl:99
 [13] macro expansion
    @ ~/ITensor Projects/FiniteT/ED_localSz_Sz_B_T.jl:133 [inlined]
 [14] macro expansion
    @ ./timing.jl:279 [inlined]
 [15] top-level scope
    @ ~/ITensor Projects/FiniteT/ED_localSz_Sz_B_T.jl:118
 [16] include(mod::Module, _path::String)
    @ Base ./Base.jl:495
 [17] exec_options(opts::Base.JLOptions)
    @ Base ./client.jl:318
 [18] _start()
    @ Base ./client.jl:552

Please know that I have 64 gig ram. I activated the term ITensors.disable_warn_order(), but I got the same error.

I’ve also experimented 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 this issue and improve the accuracy of the magnetization results in the ED method, please let me know.

4

ITensors.disable_warn_order() should disable that warning, did you execute it before running your code?

1 Like
5

Yes I added ITensors.disable_warn_order() at the beginning of the code.
But I still get memory error:

ERROR: LoadError: OutOfMemoryError()
Stacktrace:
[1] _growend!
  @ .\array.jl:1011 [inlined]
[2] push!
  @ .\array.jl:1058 [inlined]
[3] contract_blockoffsets(
    #unused#::NDTensors.AlgorithmSelection.Algorithm{:sequential, NamedTuple{(), Tuple{}}},
    boffs1::Dictionaries.Dictionary{Block{23}, Int64},
    inds1::NTuple{23, Index{Vector{Pair{QN, Int64}}}},
    labels1::NTuple{23, Int64},
    boffs2::Dictionaries.Dictionary{Block{3}, Int64},
    inds2::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}},
    labels2::Tuple{Int64, Int64, Int64},
    indsR::NTuple{24, Index{Vector{Pair{QN, Int64}}}},
    labelsR::NTuple{24, Int64})
  @ NDTensors C:\Users\PC\.julia\packages\NDTensors\bK5Ns\src\blocksparse\contract_sequential.jl:32
[4] contract_blockoffsets
  @ C:\Users\PC\.julia\packages\NDTensors\bK5Ns\src\blocksparse\contract.jl:52 [
.
.
.

I have 64 gig RAM. Please know that when I use ALPS and consider conserved total Sz, The code work without memory error.
Please also know that I edited the code as


# Calculate free energy and magnetization
function calculate_free_energy_and_magnetization(N, J1, J2xy, J2z, B_min, B_max, uB, linkdim, T)
    sites = siteinds("S=1", N; conserve_qns=true, conserve_sz=true)

    sB = Int((B_max - B_min) / uB)
    h = zeros(sB + 1)
    F_values = zeros(sB + 1)
    M_values = zeros(sB + 1)

    for (m_B, Bz) in enumerate(B_min:uB:B_max)
       h[m_B] = Bz
        statei = ["Emp" for n in 1:N]
        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_DiamondS1(N, J1, J2xy, J2z, Bz, sites)
        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

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

    return F_values, M_values, h
end

and I set parameters: N = 12 and linkdim=20.
Since, I conserved qns and sz, so the code requires statei, hence linkdim can be increased up to N^2.

6

Hi Hamid,
It sounds like you are now asking a different question from the ones above, this time about using ITensor to do ED calculations. It’s not too surprising that you are running out of memory when doing ED, which can be a memory intensive method (64 Gb is actually not very much for many-body quantum calculations).

While ITensor can be used to do simple ED calculations for benchmarking or rather small systems, it’s not really designed for state-of-the-art ED, so there will be a size where you will run out of resources quickly and I’m not surprised it’s a smaller size than ALPS, which uses specialized approaches to optimize memory usage and things in ED. ITensor is more suited for tensor network calculations like DMRG.

Best regards,
Miles

7

Dear miles,

Thanks for reply.
I successfully resolved the issue. Specifically, when employing the parameters conserve_qns=true, conserve_sz=true in ITensor, it necessitates the corresponding statei where SzT is conserved. If there exists an alternative approach wherein ITensor can automatically consider all states of the system with conserved Sz and QNs, kindly inform me.

I refined the code and managed to replicate the temperature-dependent magnetization of a Spin-1 chain using the thermal partition function Z1 += ((exp.(-E0 / T)). To achieve this, I manually inputted all states of the system with conserved Sz and QNs, and the code seems work properly.

Subsequently, I conducted a comparison between DMRG results and the ED method for a spin-1 chain with N=6, primarily to assess convergence. Despite the small chain size (N=6), the DMRG outcomes failed to converge with those of ED (see figure below).

Mag_S1_DMRG_ED

Nevertheless, I increased the number of sweeps and adjusted other DMRG parameters as follows:

linkdim = 30
sweeps = Sweeps(20)
maxdim!(sweeps, 5, 10, 20, 50, 100, 200, 300, 400, 500, 600, 800)
cutoff!(sweeps,1E-6,1E-7,1E-8,1E-9,1E-11,1E-14,1E-16,1E-18)
noise!(sweeps,1E-6,1E-7,1E-8,1E-9,1E-11,1E-14,1E-16,1E-18,0)

The discrepancy from ED persists.

I attempted to include all states with conserved Sz and QNs to encompass all energies contributing to the partition function. However, I am uncertain whether all states are being accounted for.
For instance, with the following code:

# ==========================================================
        for i in sz_min:sz_max
            #
            statei = ["Up" for n in 1:N]
            for j in 1:i
                statei[j] = "Z0"
            end
            ψi = randomMPS(sites, statei, linkdim)
            H = H_DiamondS1(N, J1, J2xy, J2z, Bz, sites)
            E0, ψ0 = dmrg(H, ψi, sweeps; write_when_maxdim_exceeds=1000)

            En = [E0]
            ψn = [ψ0]
            for i in 1:En_exited
                E, ψ = dmrg(H, ψn, ψi, sweeps; write_when_maxdim_exceeds=1000)
                println()
                push!(En, E)
                push!(ψn, ψ)
            end

            # ------------------------------------
            push!(current_energies, En)
        end
        # ==========================================================

I endeavored to encompass all states within the ‘Up’ and ‘Z0’ sectors, among others. Furthermore, I included excited states up to En_exited = 20. Despite these efforts, a noticeable discrepancy between DMRG and ED persists. Any suggestions would be greatly appreciated.
Please also know that I considered eigensolve_krylovdim=20, but it did not help to reach convergence.

Due to brevity, I omitted the portion of the code where I manually entered all states with conserved Sz and QNs. If you require it for reference, I can provide it accordingly.

8

Dear miles,

To solve the above problem concerning about the discrepancy between DMRG and ED, I just considered conserve_sz=false, and the DMRG nicely converge to the ED (see below fig).

Mag_S1_DMRG_cSzTfalse

Since the Hamiltonian commutes with Sz total, the term conserve_sz=true should work as well. I have implemented the term conserve_sz=true in another code and it works when the magnetic field is absence.
Could you please also explain that why the DMRG works well when considering conserve_sz=false in the presence of magnetic field?