Customized Operator Error: "Fluxes not all equal"

Hello ITensors Community,

I am working with the DMRG method with conserved quantum number and customized operator. Here is the simplest code for demonstration. The code exits with the error Fluxes not all equal. It works after removing the the off-diagonal term, but I indeed need to define a complicated operator.

using ITensors, ITensorMPS
ITensors.op(::OpName"T", ::SiteType"Qudit", d::Int) = 
    [1 0 1; 
    0 0 0; 
    1 0 1]
x = siteind("Qudit"; dim = 3, conserve_qns = true)  
@show op("T", x)

I sincerely appreciate it if someone could provide some insights on this problem. Thank you!

Hi @HaoyuanShi,

what exactly is your expected behavior?
I’m not sure about this, but I think the idea of conserving quantum numbers is to have operators with specific selection rules, so that when starting from a state with a given flux you can control how this is going to change. Your operator instead can either conserve the flux or increase (or decrease) it by 2.

You can bypass the error by changing the code as follows:

using ITensors, ITensorMPS

const A = [1 0 1; 0 0 0; 1 0 1]

x = siteind("S=1"; conserve_qns = true)  
T = itensor(A, x', dag(x); checkflux = false)

But, as I told you, I’m not sure if it would work later on.

Your operator conserves the parity instead, so you could try to use that as quantum number.

Example for parity qn

I’m not sure this is the proper way but it seems to work

using ITensors, ITensorMPS

alias(::SiteType"ParityQudit") = SiteType"Qudit"()

function ITensors.space(::SiteType"ParityQudit"; dim=2)
    return [QN("P", isodd(n) ? 1 : 0, 2) => 1 for n in 1:dim]
end

function ITensors.op(on::OpName, st::SiteType"ParityQudit", ds::Int...; kwargs...)
    return op(on, alias(st), ds...; kwargs...)
end

function ITensors.op(on::OpName, st::SiteType"ParityQudit", s1::Index, s_tail::Index...; kwargs...)
    rs = reverse((s1, s_tail...))
    ds = dim.(rs)
    opmat = op(on, st, ds...; kwargs...)
    return itensor(opmat, prime.(rs)..., dag.(rs)...)
end

function ITensors.op(::OpName"T", ::SiteType"ParityQudit", d::Int...; kwargs...)
    return [1 0 1;
        0 0 0;
        1 0 1]
end

x = siteind("ParityQudit"; dim=3)
T = op("T", x)

psi = MPS([1, 0, 1], (x,))
println(psi[1])
psi = apply([T], psi)
println(psi[1])

psi = MPS([0, 1, 0], (x,))
println(psi[1])
psi = apply([T], psi)
println(psi[1])

Also, if you are using spin 1 particles, consider using the S=1 site type directly.

Hello @VinceNeede

Thank you for your kind reply and detailed explanation!
Actually I am exploring the light-matter coupling with Peierls substitution, which includes the term e^{\frac {g} {\sqrt L}(a+a^\dagger)}c_j c^\dagger_{j+1} + h.c.. The exponential term is a dense matrix, and the system conserves the fermion particle number.

I encountered difficulties with DMRG due to quantum number issues at the bosonic site. Thank you for pointing out the essence of flux number. The solution is to define a new site type that has only one possible quantum number.