Liouvillian Gap in ITensor Julia

Dear ITensor developers,
I’m a beginner with ITensor Julia, and I’m trying to calculate the Liouvillian Gap for my system. However, I’ve encountered an issue with the lowest energy obtained from the non-Hermitian DMRG. It seems to be incorrect. Could someone please help me identify possible errors in my code?
Here I am using only the Ising Hamiltonian; -i(H⊗I-I⊗H^T) without any collapse operators.

using ITensors
import ITensors: op
function ITensors.space(::SiteType"Lindblad")
return 4
end

ITensors.state(::StateName"UpUp", ::SiteType"Lindblad") = [1.0, 0.0, 0.0, 0.0]
ITensors.state(::StateName"UpDn", ::SiteType"Lindblad") = [0.0, 1.0, 0.0, 0.0]
ITensors.state(::StateName"DnUp", ::SiteType"Lindblad") = [0.0, 0.0, 1.0, 0.0]
ITensors.state(::StateName"DnDn", ::SiteType"Lindblad") = [0.0, 0.0, 0.0, 1.0]

op(::OpName"SzL", ::SiteType"Lindblad") = [
0.5 0.0 0.0 0.0
0.0 -0.5 0.0 0.0
0.0 0.0 0.5 0.0
0.0 0.0 0.0 -0.5
]

op(::OpName"SzR", ::SiteType"Lindblad") = [
0.5 0.0 0.0 0.0
0.0 0.5 0.0 0.0
0.0 0.0 -0.5 0.0
0.0 0.0 0.0 -0.5
]

op(::OpName"SxL", ::SiteType"Lindblad") = [
0.0 0.5 0.0 0.0
0.5 0.0 0.0 0.0
0.0 0.0 0.0 0.5
0.0 0.0 0.5 0.0
]
op(::OpName"SxR", ::SiteType"Lindblad") = [
0.0 0.0 0.5 0.0
0.0 0.0 0.0 0.5
0.5 0.0 0.0 0.0
0.0 0.5 0.0 0.0
]

op(::OpName"SmL", ::SiteType"Lindblad") = [
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
1.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0
]

op(::OpName"SmR", ::SiteType"Lindblad") = [
0.0 0.0 0.0 0.0
1.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0
]

op(::OpName"SpL", ::SiteType"Lindblad") = [
0.0 0.0 1.0 0.0
0.0 0.0 0.0 1.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
]

op(::OpName"SpR", ::SiteType"Lindblad") = [
0.0 1.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 1.0
0.0 0.0 0.0 0.0
]
N = 8
sites = siteinds(“Lindblad”,N)
opsum = OpSum()

add!(opsum,-4im,“SzL”,1,“SzL”,2)

opsum += (4im,“SzR”,1,“SzR”,2)
opsum += (-4im,“SzL”,2,“SzL”,3)
opsum += (4im,“SzR”,2,“SzR”,3)
opsum += (-4im,“SzL”,3,“SzL”,4)
opsum += (4im,“SzR”,3,“SzR”,4)
H = MPO(opsum,sites)

psi0 = randomMPS(sites,10)
nsweeps = 5
maxdim = [20,40,80,100,120]
cutoff = [1E-10]

energy, psi = dmrg(H,psi0; nsweeps, maxdim, cutoff, ishermitian=false)

@show energy[0]

1 Like