entanglement entropy in tdvp calculation

Dear iTensorJulia Community,

I am encountering an issue with the attached code snippet, which I borrowed from the TDVP example and extended with a function to calculate entanglement entropy. However, upon running the code, I am encountering an error, as shown in the attached files.

function main()
  
    function entropy_von_neumann(psi::MPS, b::Int)
      s = siteinds(psi)  
      orthogonalize!(psi, b)
      _,S = svd(psi[b], (linkind(psi, b-1), s[b]))
      SvN = 0.0
      for n in 1:dim(S, 1)
        p = S[n,n]^2
        SvN -= p * log(p)
      end
      return SvN
    end
  
  
    function heisenberg(N)
      os = OpSum()
      for j in 1:(N - 1)
        os += 0.5, "S+", j, "S-", j + 1
        os += 0.5, "S-", j, "S+", j + 1
        os += "Sz", j, "Sz", j + 1
      end
      return os
    end
    
    N = 10
    cutoff = 1e-12
    tau = 0.1
    ttotal = 5.0
  
    s = siteinds("S=1/2", N; conserve_qns=true)
    H = MPO(heisenberg(N), s)
  
    function step(; sweep, bond, half_sweep)
      if bond == 1 && half_sweep == 2
        return sweep
      end
      return nothing
    end
  
    function current_time(; current_time, bond, half_sweep)
      if bond == 1 && half_sweep == 2
        return current_time
      end
      return nothing
    end
  
    function measure_sz(; psi, bond, half_sweep)
      if bond == 1 && half_sweep == 2
        return expect(psi, "Sz"; sites=1)
      end
      return nothing
    end
  
  
    function measure_SvN(; psi, bond, half_sweep)
      if bond == 1 && half_sweep == 2
        return entropy_von_neumann(psi, 1)
      end
      return nothing
    end
  
  
    function return_state(; psi, bond, half_sweep)
      if bond == 1 && half_sweep == 2
        return psi
      end
      return nothing
    end
    
  
    obs = observer(
      "steps" => step, "times" => current_time, 
      "psis" => return_state, "Sz" => measure_sz,
      "SvN" => measure_SvN,
    )
  
    psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")
  
    @show psi
  
    
    psi_f = tdvp(
      H,
      -im * ttotal,
      psi;
      time_step=-im * tau,
      cutoff,
      outputlevel=3,
      normalize=false,
      (observer!)=obs,
    )
  
    @show psi_f
    
    
    steps = obs.steps
    times = obs.times
    psis  = obs.psis
    Sz    = obs.Sz
    SvN    = obs.SvN
   
    
  
    println("\nResults")
    println("=======")
    for n in 1:length(steps)
      print("step = ", steps[n])
      print(", time = ", round(times[n]; digits=3))
      print(", |⟨ψⁿ|ψⁱ⟩| = ", round(abs(inner(psis[n], psi));   digits=3))
      print(", |⟨ψⁿ|ψᶠ⟩| = ", round(abs(inner(psis[n], psi_f)); digits=3))
      print(", ⟨Sᶻ⟩ = ", round(Sz[n]; digits=3))
      print(", SvN = ",  round(entropy_von_neumann(psis[n], 1); digits=3))
      println()
    end
  
  
    magz = expect(psi_f,"Sz")
    for (j,mz) in enumerate(magz)
        println("$j $mz")
    end
  

  
    return nothing
  end
  
  main()

here is the ERROR:
###################
###################

MethodError: no method matching _indices(::Tuple{}, ::Nothing)

Closest candidates are: _indices(::Tuple, !Matched::Vector) @ ITensors ~/.julia/packages/ITensors/WMeVS/src/indexset.jl:55 _indices(!Matched::Vector, ::Any) @ ITensors ~/.julia/packages/ITensors/WMeVS/src/indexset.jl:52 _indices(::Any, !Matched::Vector) @ ITensors ~/.julia/packages/ITensors/WMeVS/src/indexset.jl:53 … Stacktrace: [1] (::Base.BottomRF{typeof(ITensors._indices)})(acc::Tuple{}, x::Nothing) @ Base ./reduce.jl:86 [2] afoldl(::Base.BottomRF{typeof(ITensors._indices)}, ::Tuple{}, ::Nothing, ::Index{Vector{Pair{QN, Int64}}}) @ Base ./operators.jl:544 [3] _foldl_impl(op::Base.BottomRF{typeof(ITensors._indices)}, init::Tuple{}, itr::Tuple{Nothing, Index{Vector{Pair{QN, Int64}}}}) @ Base ./reduce.jl:68 [4] foldl_impl(op::Base.BottomRF{typeof(ITensors._indices)}, nt::Tuple{}, itr::Tuple{Nothing, Index{Vector{Pair{QN, Int64}}}}) @ Base ./reduce.jl:48 [5] mapfoldl_impl(f::typeof(identity), op::typeof(ITensors._indices), nt::Tuple{}, itr::Tuple{Nothing, Index{Vector{Pair{QN, Int64}}}}) @ Base ./reduce.jl:44 [6] mapfoldl(f::Function, op::Function, itr::Tuple{Nothing, Index{Vector{Pair{QN, Int64}}}}; init::Tuple{}) @ Base ./reduce.jl:175

@ ~/.julia/packages/ITensorTDVP/njU8f/src/tdvp.jl:47 [inlined] [39] main() @ Main ~/Desktop/FIles/mehdiAbdi/Ising-longRange/test-model/tdmrg-itensor-Julia/main.ipynb:83 [40] top-level scope @ ~/Desktop/FIles/mehdiAbdi/Ising-longRange/test-model/tdmrg-itensor-Julia/main.ipynb:128

###################
###################

From my observations, it appears that the output of TDVP is not in the MPS format, which is causing the error. Strangely though, when I use the TDVP output and pass it to calculate the “Sz” expectation, as seen towards the end of the code, it works fine. It seems that the built-in function “expect()” expects the input to be in the MPS format, but my function for entanglement entropy isn’t working with it.

Any assistance regarding this issue would be greatly appreciated.

Best regards,
Javad

Try changing the line:

      _,S = svd(psi[b], (linkind(psi, b-1), s[b]))

to:

      _,S = svd(psi[b], (linkinds(psi, b-1)..., s[b]))

Basically what is happening is that for b == 1, you are trying to grab linkind(psi, 0), which doesn’t exist so it returns nothing. linkinds(psi, b-1) outputs a collection of link indices, so is either an empty collection ([]) if b == 1 or a collection of length one containing the link index ([linkind(psi, b - 1)]) if b > 1 (assuming there is only one link index, there could in principle be more than one link index though that is not standard and you would have to construct that yourself). Changing that line to (linkinds(psi, b-1)..., s[b]) means that if b == 1 it makes a collection of indices just containing the site index, (s[b],), while if b > 1 it makes the collection equal to (linkind(psi, b-1), s[b]). For example:

julia> using ITensors

julia> s = siteinds("S=1/2", 4);

julia> psi = MPS(s, n -> isodd(n) ? "Up" : "Dn");

julia> linkind(psi, 0)

julia> (linkind(psi, 0), s[1])
(nothing, (dim=2|id=100|"S=1/2,Site,n=1"))

julia> (linkinds(psi, 0)..., s[1])
((dim=2|id=100|"S=1/2,Site,n=1"),)

julia> linkind(psi, 1)
(dim=1|id=471|"Link,l=1")

julia> (linkind(psi, 1), s[2])
((dim=1|id=471|"Link,l=1"), (dim=2|id=541|"S=1/2,Site,n=2"))

julia> (linkinds(psi, 1)..., s[2])
((dim=1|id=471|"Link,l=1"), (dim=2|id=541|"S=1/2,Site,n=2"))

The svd code doesn’t understand when you input the indices (nothing, (dim=2|id=100|"S=1/2,Site,n=1")) in the b == 1 case so it is throwing that opaque error.

I see this issue is coming up a lot (Entanglement Entropy Bug, About measuring entropy, etc.). I think everyone is copying the code from the docs here: MPS and MPO Examples · ITensors.jl which should be updated to the way I suggested. I’ll update that now.

1 Like

Dear mtfishman,

great and many thanks for your reply.