Inconsistent behaviour of combiner with and without qns

ITensor Julia Questions
1

I tried constructing a combiner for an index and its dagger to obtain the doubled space of a fermion (where a vectorized density matrix would live). Depending on whether I use quantum numbers or not, I get different results. I think this is because the combiner has to reorder in the presence of quantum numbers, but not in their absence.

For the specific example below, I can define a workaround by enforcing the reorder (see last line), but generally this looks dangerous and unexpected. The example below is not minimal, but gives a bit more context


using ITensors
function ITensors.space(::SiteType"FoldedFermion"; conserve_qns::Bool=false)
    conserve_qns && return [QN("Nf",-1,-1)=>1,QN("Nf",0,-1)=>2,QN("Nf",+1,-1)=>1]
    return 4
end
function folder(;conserve_qns::Bool=false)
    f = ITensors.siteind("Fermion"; conserve_qns)
    cmb = combiner(f,dag(f)')
    ff = ITensors.siteind("FoldedFermion"; conserve_qns)
    cmb *= denseblocks(delta(ff ,dag(combinedind(cmb)))) # replace index and make dense
    return cmb
end

@show dense(folder(;conserve_qns=true)) # correct
@show fldr = folder(;conserve_qns=false) # wrong: different index order from conserved_qns
@show fldr = permute(fldr,inds(fldr)[[2,3,1]]) # wrong: same index order but different results from conserved_qns
f = ind(fldr,1)
@show permute(prime(ITensor([0 1; 1 0],f', f'') * fldr, -1; plev=2), inds(fldr)) # correct: apply X to outgoing (primed) index

Then I would like to created the (vectorized) density matrix of an occupied fermionic site and convert it to an actual density matrix (or vice-versa), but definitions become inconsistent in the non-conserving case

ITensors.state(::StateName"1",::SiteType"FoldedFermion") = [0,0,1,0]

# occupied fermion density matrix: outputs [0 0; 0 1]
fldr = folder(;conserve_qns=true)
@show dense(ITensors.state(inds(fldr; tags="FoldedFermion")[end],"1") * dag(fldr)) 

# wrong: outputs [0 1; 0 0]
fldr = folder(;conserve_qns=false)
@show dense(ITensors.state(inds(fldr; tags="FoldedFermion")[end],"1") * dag(fldr))

Is there a way to fix the output basis order of a combiner? Or maybe you have some better suggestion?

Thanks,
Matan