# Problem using TEBD algorithm to represent Gibbs states

**URL:** <https://itensor.discourse.group/t/problem-using-tebd-algorithm-to-represent-gibbs-states/2356>\
**Category:** General Discussion\
**Tags:** julia, mpo\
**Created:** [April 23, 2025, 2:17pm UTC](https://itensor.discourse.group/t/problem-using-tebd-algorithm-to-represent-gibbs-states/2356 "2025-04-23T14:17:30Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![jorgesanseg](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/jorgesanseg/32/684_2.png) [@jorgesanseg](https://itensor.discourse.group/u/jorgesanseg)\
**Post date:** [April 23, 2025, 2:17pm UTC](https://itensor.discourse.group/t/problem-using-tebd-algorithm-to-represent-gibbs-states/2356/1 "2025-04-23T14:17:30Z")

</div>

Hi,

I am trying to get thermal states \rho\_\beta = e^{-\beta H} mean values of local observables for long-range interacting system. More concretely the Long Range Spin chain:

H\_{\text{LR, Ising}}= J \sum\_{i \neq j}^N \frac{\sigma\_i^x \sigma\_{j}^x}{|i-j|^{\alpha}} + h \sum\_{i=1}^N \sigma\_i^z

Notice that the smaller the \alpha, the “least” local the model is, and thereby MPs and MPOs representations are less efficient as the entanglement of the system is greater.

We construct the time evolution operator of the LR-Ising model explicitly in the short-range case([Radware Bot Manager Captcha](http://dx.doi.org/10.1088/1367-2630/12/2/025012)) and then applying a super-extensive amount of swap gates to exactly simulate the power-law interactions in a finite system size. This yields a Trotterized representation of the time evolution operator U(\delta t) = \mathrm{e}^{-iH \delta t} in form of an MPO. We concatenate the Trotterized time evoltion to yield the desired final _imaginary time_ or \beta U(t) \simeq \prod\_{i=1}^{n\_t} U(\delta t), with t = n\_t \cdot \delta t.

When I perform the simulations for \alpha \leq 1.5, despite the bond dimension not to collapse the evolution stops and we have the following error:

ERROR: LoadError: BoundsError: attempt to access 0-element Vector{Pair{QN, Int64}} at index [1]  
Stacktrace:  
[1] throw\_boundserror(A::Vector{Pair{QN, Int64}}, I::Tuple{Int64})  
@ Base ./essentials.jl:14  
[2] getindex  
@ ./essentials.jl:916 [inlined]  
[3] combineblocks(qns::Vector{Pair{QN, Int64}})  
@ ITensors ~/.julia/packages/ITensors/4M2ep/src/qn/qnindex.jl:413  
[4] combineblocks  
@ ~/.julia/packages/ITensors/4M2ep/src/qn/qnindex.jl:475 [inlined]  
[5] combiner(inds::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}; dir::Nothing, tags::String)  
@ ITensors ~/.julia/packages/ITensors/4M2ep/src/qn/qnitensor.jl:444  
[6] combiner(::Index{Vector{Pair{QN, Int64}}}, ::Vararg{Index{Vector{Pair{QN, Int64}}}}; kwargs::@Kwargs{dir::Nothing})  
@ ITensors ~/.julia/packages/ITensors/4M2ep/src/tensor\_operations/itensor\_combiner.jl:7  
[7] svd(A::ITensor, Linds::Vector{Index{Vector{Pair{QN, Int64}}}}; leftdir::Nothing, rightdir::Nothing, lefttags::TagSet, righttags::TagSet, mindim::Int64, maxdim::Int64, cutoff::Float64, alg::Nothing, use\_absolute\_cutoff::Nothing, use\_relative\_cutoff::Nothing, min\_blockdim::Nothing, utags::TagSet, vtags::TagSet)  
@ ITensors ~/.julia/packages/ITensors/4M2ep/src/tensor\_operations/matrix\_decomposition.jl:153  
[8] svd  
@ ~/.julia/packages/ITensors/4M2ep/src/tensor\_operations/matrix\_decomposition.jl:110 [inlined]  
[9] factorize\_svd(A::ITensor, Linds::Vector{Index{Vector{Pair{QN, Int64}}}}; singular\_values!::Nothing, ortho::String, alg::Nothing, dir::Nothing, mindim::Int64, maxdim::Int64, cutoff::Float64, tags::TagSet, use\_absolute\_cutoff::Nothing, use\_relative\_cutoff::Nothing, min\_blockdim::Nothing)  
@ ITensors ~/.julia/packages/ITensors/4M2ep/src/tensor\_operations/matrix\_decomposition.jl:614  
[10] factorize(A::ITensor, Linds::Vector{Index{Vector{Pair{QN, Int64}}}}; mindim::Int64, maxdim::Int64, cutoff::Float64, ortho::String, tags::TagSet, plev::Nothing, which\_decomp::Nothing, eigen\_perturbation::Nothing, svd\_alg::Nothing, use\_absolute\_cutoff::Nothing, use\_relative\_cutoff::Nothing, min\_blockdim::Nothing, singular\_values!::Nothing, dir::Nothing)  
@ ITensors ~/.julia/packages/ITensors/4M2ep/src/tensor\_operations/matrix\_decomposition.jl:808  
[11] contract(::NDTensors.BackendSelection.Algorithm{:zipup, @NamedTuple{}}, A::MPO, B::MPO; cutoff::Float64, maxdim::Int64, mindim::Int64, kwargs::@Kwargs{})  
@ ITensors.ITensorMPS ~/.julia/packages/ITensors/4M2ep/src/ITensorMPS/mpo.jl:845  
[12] contract  
@ ~/.julia/packages/ITensors/4M2ep/src/ITensorMPS/mpo.jl:815 [inlined]  
[13] contract(A::MPO, B::MPO; alg::String, kwargs::@Kwargs{cutoff::Float64, maxdim::Int64})  
@ ITensors.ITensorMPS ~/.julia/packages/ITensors/4M2ep/src/ITensorMPS/mpo.jl:808  
[14] contract  
@ ~/.julia/packages/ITensors/4M2ep/src/ITensorMPS/mpo.jl:807 [inlined]  
[15] product(A::MPO, B::MPO; kwargs::@Kwargs{cutoff::Float64, maxdim::Int64})  
@ ITensors.ITensorMPS ~/.julia/packages/ITensors/4M2ep/src/ITensorMPS/mpo.jl:885  
[16] top-level scope  
@ ~/jssegovia/Codes/CorrelationOfDecayLR/CorrelationOfDecayLR/GibbsExpectationValue\_tosend.jl:87  
in expression starting at /lustre/home/ift/jssegovia/jssegovia/Codes/CorrelationOfDecayLR/CorrelationOfDecayLR/GibbsExpectationValue\_tosend.jl:83

The code is the following, and I’ve used the parameters: \delta t=0.01, SVDcutoff=1e-18, and the evolution stopped at \beta=2.85 when I was aiming \beta=3. The problem gets worse when decreasing \alpha:

The code is something like this:

```auto
for (i, t) in enumerate(β_values)
    contraction_U[:] = apply(U_δt, contraction_U; cutoff=SVDcutoff, maxdim=SVDmaxdim) #here U_δt is the infinitesimal imaginary time evolution MPO, and contraction_U is initialized as the Identity MPO

    trace_values[i] = tr(contraction_U)
    A_values[i]= tr(apply( contraction_U, A; cutoff=SVDcutoff, maxdim=SVDmaxdim))*(1/trace_values[i])

    maxdim = maxlinkdim(contraction_U)

    open(csv_filename, "a") do file
        writedlm(file, [t trace_values[i] A_values[i] maxdim], ',')
    end

    @printf("After t=%g\tmaxlinkdim=%d\ttimestep %d out of %d", t, maxdim, i, ntimes)
    println(" runtime: $(Dates.canonicalize(Dates.now() - t1))")
end

```

I don’t have a good intuition of what is going on, If someone has it would be really helpful!

---

<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:** [April 23, 2025, 4:22pm UTC](https://itensor.discourse.group/t/problem-using-tebd-algorithm-to-represent-gibbs-states/2356/2 "2025-04-23T16:22:25Z")

</div>

It would be difficult to debug this without a fully runnable code.

Can you share an example of MPOs you are producting together that leads to that error (say by saving them with HDF5 or JLD2) along with the values of `cutoff` and `maxdim` you passed to `apply`?

---

<div class="post-metadata">

**Author:** ![jorgesanseg](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/jorgesanseg/32/684_2.png) [@jorgesanseg](https://itensor.discourse.group/u/jorgesanseg)\
**Post date:** [April 23, 2025, 4:43pm UTC](https://itensor.discourse.group/t/problem-using-tebd-algorithm-to-represent-gibbs-states/2356/3 "2025-04-23T16:43:01Z")

</div>

Sure! The functions I use for the MPO representation are :

```auto
function build_expH_ising_murg(sites, JXX::Real, hz::Real, dt::Number)

  

  # For real dt this does REAL time evolution
  # I should have already taken into account both the - sign in exp(-iHt)
  # and the overall minus in Ising H= -(JXX+hZ)

  # dt = JXX * dt 
  # JXX= JXX * dt ## why?

  cosg = cos(hz * dt * 0.5)
  sing = sin(hz * dt * 0.5)

  N = length(sites)
  U_t = MPO(N)

  link_dimension = 2

  # linkindices = hasqns(sites) ? [Index(link_dimension, "Link,l=$(n-1)") for n in 1:N+1] : [Index(link_dimension, "Link,l=$(n-1)") for n in 1:N+1]
  linkindices = [Index(link_dimension, "Link,l=$(n-1)") for n in 1:N+1]

  for n in 1:N
    # siteindex s
    s = sites[n]

    # left link index ll with daggered QN conserving direction (if applicable)
    ll = dag(linkindices[n])
    # right link index rl
    rl = linkindices[n+1]

    combz = (1 - 2 * sing^2) * op(sites, "Id", n) + im * 2 * sing * cosg * op(sites, "Z", n)
    X = op(sites, "X", n)

    if n == 1
      #U_t[n] = ITensor(ComplexF64, dag(s), s', dag(rl))
      U_t[n] = onehot(rl => 1) * sqrt(cos(JXX*dt)) * combz
      U_t[n] += onehot(rl => 2) * sqrt(im * sin(JXX*dt)) * X
    elseif n == N
      #U_t[n] = ITensor(ComplexF64, ll, dag(s), s')
      U_t[n] = onehot(ll => 1) * sqrt(cos(JXX*dt)) * combz
      U_t[n] += onehot(ll => 2) * sqrt(im * sin(JXX*dt)) * X

    else
      #U_t[n] = ITensor(ComplexF64, ll, dag(s), s', dag(rl))

      U_t[n] = onehot(ll => 1, rl => 1) * cos(JXX*dt) * combz
      U_t[n] += onehot(ll => 1, rl => 2) * sqrt(im * sin(JXX*dt)) * sqrt(cos(JXX*dt)) * X
      U_t[n] += onehot(ll => 2, rl => 1) * sqrt(im * sin(JXX*dt)) * sqrt(cos(JXX*dt)) * X
      U_t[n] += onehot(ll => 2, rl => 2) * im * sin(JXX*dt) * combz
    end
  end
  return U_t
end

function XX_NN_layer(sites, J::Float64)
  N = length(sites)
  return [op("Rxx", sites[j], sites[j+1], ϕ=-J) for j in 1:N-1 ]
end

function XX_powerlaw_layer(sites, α::Float64, J::Number)
  N = length(sites)
  layer = Vector{ITensor}(undef,0)
  # number_of_terms = 3 * N * (N - 1)/2 + 3*(N-1)
  for i in 1:N-1
    for j in i+1:N-1
      # verify factor!
      push!(layer, op("Rxx", sites[j - 1], sites[j], ϕ=-J * (1.0 / ((j - i)^α))) )
      push!(layer, op("Swap", sites[j - 1], sites[j]))
    end
    push!(layer, op("Rxx", sites[N - 1], sites[N], ϕ=-J * (1.0 / ((N - i)^α))))
    push!(layer, [op("Swap", sites[j - 1], sites[j]) for j in N-1:-1:i+1]...)
  end
  return layer
end

function XX_NTNT_layer(sites, J::Float64)
  N = length(sites)
  layer = Vector{ITensor}(undef,0)
  # number_of_terms = 3 * N * (N - 1)/2 + 3*(N-1)
  for i in 1:N-2
    push!(layer, op("Rxx", sites[i], sites[i+1], ϕ=-J))
    push!(layer, op("Swap", sites[i], sites[i+1]))
    push!(layer, op("Rxx", sites[i+1], sites[i+2], ϕ=-J))
    push!(layer, op("Swap", sites[i+1], sites[i]))
  end
  push!(layer, op("Rxx", sites[N-1], sites[N], ϕ=-J))
  return layer
end

function Rz_rotation_layer(sites, ϕ::Number)
  N = length(sites)
  # verify the factor!!!
  return MPO([op("Rz", sites[j], θ=-2ϕ) for j in 1:N], 0, N+1)
end

function Rx_rotation_layer(sites, ϕ::Number)
  N = length(sites)
  # verify the factor!!!
  return MPO([op("Rx", sites[j], θ=-2ϕ) for j in 1:N], 0, N+1)
end

function timeEvo_Cirac_LR_TFI(sites, α::Float64, J::Float64,hp::Float64,δtn::Number)
  Uz = Rz_rotation_layer(sites,hp*δtn*0.5)
  Uxx = XX_powerlaw_layer(sites,α,J*δtn)

  Uz_Uxx = mapreduce(
    e -> apply(e, Uz; maxdim=SVDmaxdim,cutoff=SVDcutoff),
    (previousResult, newElement) -> apply(newElement, previousResult; maxdim=SVDmaxdim, cutoff=SVDcutoff),
  Uxx
  )

  return mapreduce(e -> apply(e, Uz_Uxx; maxdim=SVDmaxdim,cutoff=SVDcutoff),
  (previousResult, newElement) -> apply(newElement, previousResult; maxdim=SVDmaxdim,cutoff=SVDcutoff),
  Uz)
end

function timeEvo_Cirac_LR_TFI_corr(sites, α::Float64, J::Float64,hp::Float64,δtn::Number)
  Uz = Rz_rotation_layer(sites,hp*δtn*0.5)
  Uxx = XX_powerlaw_layer(sites,α,J*δtn)
  return mapreduce(e -> apply(Uz,e; maxdim=SVDmaxdim,cutoff=SVDcutoff), e-> apply(Uxx,e; maxdim=SVDmaxdim,cutoff=SVDcutoff), e-> apply(Uz,e; maxdim=SVDmaxdim,cutoff=SVDcutoff))

  # Uz_Uxx = mapreduce(
  # e -> apply(Uz,e; maxdim=SVDmaxdim,cutoff=SVDcutoff),
  # (previousResult, newElement) -> apply(newElement, previousResult; maxdim=SVDmaxdim, cutoff=SVDcutoff),
  # Uxx
  # )

  # return mapreduce(e -> apply(Uz_Uxx,e; maxdim=SVDmaxdim,cutoff=SVDcutoff),
  # (previousResult, newElement) -> apply(newElement, previousResult; maxdim=SVDmaxdim,cutoff=SVDcutoff),
  # Uz)
end

function timeEvo_Cirac_LR_TFI_NewVersion(sites, α::Float64, J::Float64,hp::Float64,δtn::Number,SVDcutoff::Float64,SVDmaxdim::Int64)
  OperatorMPO = Rz_rotation_layer(sites,hp*δtn*0.5)
  Uz = Rz_rotation_layer(sites,hp*δtn*0.5)
  Uxx = XX_powerlaw_layer(sites,α,J*δtn)
  
  # OperatorMPO[:] = apply(Uz, OperatorMPO; cutoff=SVDcutoff, maxdim=SVDmaxdim)
  OperatorMPO[:] = apply(Uxx, OperatorMPO; cutoff=SVDcutoff, maxdim=SVDmaxdim)
  OperatorMPO[:] = apply(Uz, OperatorMPO; cutoff=SVDcutoff, maxdim=SVDmaxdim)

  return OperatorMPO
end

```

and the operator is defined as:

```auto
const U_δt = timeEvo_Cirac_LR_TFI_NewVersion(sites, α, J, hp, δβ, SVDcutoff, SVDmaxdimOperator)
const contraction_U = MPO(map(n -> op(sites, "Id", n), 1:length(sites)))

```

Then, the parameters you asked for: const SVDmaxdimOperator = 30, const SVDcutoff = 1e-18, const SVDmaxdim=1800, const \delta \beta = - 0.01\* im, const J=1, const hp=0.5. For \alpha\leq 1.5 system size above 55 I start having this problem, and the evolution does not finish.

pd: I am creating a hdf5 file with the contraction, I did not have it before as I simply loaded the data in a csv to save some space

Thank you very much for your help!

---

<div class="post-metadata">

**Author:** ![jorgesanseg](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/jorgesanseg/32/684_2.png) [@jorgesanseg](https://itensor.discourse.group/u/jorgesanseg)\
**Post date:** [May 29, 2025, 9:45am UTC](https://itensor.discourse.group/t/problem-using-tebd-algorithm-to-represent-gibbs-states/2356/4 "2025-05-29T09:45:25Z")

</div>

Hi, sorry for the delay, but I finally have a H5 file with the MPO during imaginary time evolution. I’ve found no other way to send it than sharing this link with a zip: [https://drive.google.com/file/d/1D51dqBeY0oj\_XaoMJLxxB0FcgbFL8N5C/view?usp=drive\_link](https://drive.google.com/file/d/1D51dqBeY0oj_XaoMJLxxB0FcgbFL8N5C/view?usp=drive_link)

If I’d rather send it differently let me know, thanks for your time!