# Conserve_qns not working for custom MPO

+1 vote
edited

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

+1 vote

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)
space = hasqns(rho[1]) ? (QN()=>1) : 1
for n=1:N
end
for n=1: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