Conserve_qns not working for custom MPO

Question

Hey everyone,

I am performing Heisenberg evolution of an MPO. I have performed several weeks to make sure the code is giving the correct results based off benchmarks to other works, and now I want to make my code more efficient before adding multi-threading. I figured if I could use the conserved_qns = true feature the MPO would more sparse and speed up computation. It seems the issue when my gates(TEBD) apply to the MPO I get an error.

function initial(N::Int, mu::Float64, index)
      rho = MPO(N)
      for x=1:Int(0.5*N)
           rho[x] = 0.5*(op("Id", index[x]) + 2*mu*op("Sz",index[x])
      end
      for x=Int(0.5*N+1):N
           rho[x] = 0.5*(op("Id", index[x]) - 2*mu*op("Sz",index[x])
       end
       return rho
    end

The gate function is the same the one found in the examples on the GitHub page. All functions work properly when the conserve_qns=true, the issue occurs when using "apply" to act on the MPO.

apply(gates, rho; cutoff=1e-8, maxdim=64, apply_dag=true)

When this is done I get the error, "In 'svd' the left or right indices are empty.." I am not sure why this occurs when the siteinds are set to conserve but any help is appreciated.

Chris

Answer 1

Hi Chris,
I found the issue, which is that our apply function fails in certain cases if the MPO passed to it does not have any link or virtual indices. Ideally it would be able to just work in that case, since I'm sure you'd agree it's more convenient to just make the MPO the way you did without putting in link indices by hand also.

But for now, I've made an improved version of your "initial" function as well as a helper function called "putlinks!" which puts dimension-1 link indices into the MPO. After using putlinks!, I find that the apply code runs fine even with quantum numbers (basing it on the example code here: https://github.com/ITensor/ITensors.jl/blob/main/examples/gate_evolution/mpo_gate_evolution.jl )

function putlinks!(rho::MPO)
  N = length(rho)
  links = Index[]
  space = hasqns(rho[1]) ? (QN()=>1) : 1
  for n=1:N
    l = Index(space;tags="Link,n=$n")
    push!(links,l)
  end
  for n=1:N-1
    rho[n] *= dag(onehot(links[n]=>1))
    rho[n+1] *= onehot(links[n]=>1)
  end
end

function initial(N::Int, mu::Float64, index)
  rho = MPO(N)
  Nh = div(N,2)
  for x=1:Nh
    rho[x] = 0.5*(op("Id", index[x]) + 2*mu*op("Sz",index[x]))
  end
  for x=Nh+1:N
    rho[x] = 0.5*(op("Id", index[x]) - 2*mu*op("Sz",index[x]))
  end
  putlinks!(rho)
  return rho
end

Comments

Wow, this is a lot better than just doing it by hand to get things to work. I thought it was the apply function, I tracked the bug into that section of code. I am glad you were able to figure it out. Hopefully, it will lead to future improvements!

thanks agian