# Conserving total quantum number SzT in spinful fermionic ladder

**URL:** <https://itensor.discourse.group/t/conserving-total-quantum-number-szt-in-spinful-fermionic-ladder/1083>\
**Category:** ITensor Julia Questions\
**Tags:** dmrg\
**Created:** [August 9, 2023, 11:57am UTC](https://itensor.discourse.group/t/conserving-total-quantum-number-szt-in-spinful-fermionic-ladder/1083 "2023-08-09T11:57:10Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![hamid](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/hamid/32/236_2.png) [@hamid](https://itensor.discourse.group/u/hamid)\
**Post date:** [August 9, 2023, 11:57am UTC](https://itensor.discourse.group/t/conserving-total-quantum-number-szt-in-spinful-fermionic-ladder/1083/1 "2023-08-09T11:57:10Z")

</div>

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:

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

```

I initialize the state `statei` as:

```auto
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:

```auto
# 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:

```auto
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.

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [August 10, 2023, 1:04pm UTC](https://itensor.discourse.group/t/conserving-total-quantum-number-szt-in-spinful-fermionic-ladder/1083/2 "2023-08-10T13:04:52Z")

</div>

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](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.

---

<div class="post-metadata">

**Author:** ![hamid](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/hamid/32/236_2.png) [@hamid](https://itensor.discourse.group/u/hamid)\
**Post date:** [February 1, 2024, 1:28pm UTC](https://itensor.discourse.group/t/conserving-total-quantum-number-szt-in-spinful-fermionic-ladder/1083/3 "2024-02-01T13:28:32Z")

</div>

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:

```auto
# 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:

```auto
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.

---

<div class="post-metadata">

**Author:** ![mtfishman](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/mtfishman/32/11_2.png) [@mtfishman](https://itensor.discourse.group/u/mtfishman)\
**Post date:** [February 1, 2024, 2:22pm UTC](https://itensor.discourse.group/t/conserving-total-quantum-number-szt-in-spinful-fermionic-ladder/1083/4 "2024-02-01T14:22:32Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![hamid](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/hamid/32/236_2.png) [@hamid](https://itensor.discourse.group/u/hamid)\
**Post date:** [February 2, 2024, 8:00am UTC](https://itensor.discourse.group/t/conserving-total-quantum-number-szt-in-spinful-fermionic-ladder/1083/5 "2024-02-02T08:00:55Z")

</div>

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

```auto
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

```auto

# 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`.

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [February 3, 2024, 4:00pm UTC](https://itensor.discourse.group/t/conserving-total-quantum-number-szt-in-spinful-fermionic-ladder/1083/6 "2024-02-03T16:00:58Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![hamid](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/hamid/32/236_2.png) [@hamid](https://itensor.discourse.group/u/hamid)\
**Post date:** [February 6, 2024, 7:36pm UTC](https://itensor.discourse.group/t/conserving-total-quantum-number-szt-in-spinful-fermionic-ladder/1083/7 "2024-02-06T19:36:03Z")

</div>

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](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/5433268c2074d7756aa83c343983795bd4b8be01.png)

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

```auto
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:

```auto
# ==========================================================
        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.

---

<div class="post-metadata">

**Author:** ![hamid](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/hamid/32/236_2.png) [@hamid](https://itensor.discourse.group/u/hamid)\
**Post date:** [February 7, 2024, 9:40pm UTC](https://itensor.discourse.group/t/conserving-total-quantum-number-szt-in-spinful-fermionic-ladder/1083/8 "2024-02-07T21:40:10Z")

</div>

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](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/c1787270bd90ff6890209c3986bb4f97b4d7e3f1.png)

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?