“Eigen currently only supports block diagonal matrices.” when applying an MPO to a MPS with QNs

[phx]

I’m trying to find a way to write the mixed state of a many-body system of spinless fermions as an MPS of the Electron site type. The idea is essentially to map the primed site index of each MPS tensor to the ↑ degree of freedom of the electron, and the dag’ed site index to ↓.

To do so, I wrote the following function, which should be the “inverse” of the electron_to_fermion_site function from Best way to move from electron to fermion site type - #3 by corbett5. I need to deal with QNs since the fermion system I am studying enjoys the conservation of the particle number (that translates to the number conservation of electrons, too).

function fermion_to_electron_site(t::ITensor, fup::Index, fdn::Index, es::Index)
    nonSiteInds = [i for i in inds(t) if (i != fup && i != fdn)]
    t = permute(t, nonSiteInds..., fup, fdn)
    newInds = copy(nonSiteInds)
    push!(newInds, es)
    newT = ITensors.BlockSparseTensor(eltype(t), Block{length(newInds)}[], newInds)
    for coords in
        CartesianIndex([1 for _ in nonSiteInds]...):CartesianIndex(
        [dim(i) for i in nonSiteInds]...
    )
        newT[coords.I..., 1] = t[coords.I..., 1, 1]
        newT[coords.I..., 2] = t[coords.I..., 2, 1]
        newT[coords.I..., 3] = t[coords.I..., 1, 2]
        newT[coords.I..., 4] = t[coords.I..., 2, 2]
    end
    return itensor(newT)
end

This function takes an MPS x representing the pure state of a Fermion system and returns the corresponding density matrix, as an Electron system.

function projector_sf(x::MPS; kwargs...)
    r = projector(x)
    se = siteinds("Electron", length(x); kwargs...)
    return MPS([fermion_to_electron_site(r[k], siteinds(r,k)..., se[k]) for k in 1:length(r)])
end

Now, I’m trying to test this construction by computing the expectation value of some observables both on x and on projector_sf(x), and seeing if I get the same answer.
However, I can’t even apply an MPO to the second MPS. This code, in fact,

using ITensors, ITensorMPS
N = 8
sf = siteinds("Fermion", N; conserve_nf=true)
x = MPS(sf, "Occ")
xs = projector_sf(x; conserve_nf=true)
ops = OpSum()
ops += "n↑", 1, "n↑", 2
apply(MPO(ops, siteinds(xs)), xs)

yields the following error:

ERROR: Eigen currently only supports block diagonal matrices.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] eigen(T::LinearAlgebra.Hermitian{…}; min_blockdim::Nothing, mindim::Int64, maxdim::Int64, cutoff::Float64, use_absolute_cutoff::Nothing, use_relative_cutoff::Nothing)
    @ NDTensors ~/.julia/packages/NDTensors/6HwdB/src/blocksparse/linearalgebra.jl:230
  [3] eigen(A::ITensor, Linds::Vector{…}, Rinds::Vector{…}; mindim::Int64, maxdim::Int64, cutoff::Float64, use_absolute_cutoff::Nothing, use_relative_cutoff::Nothing, ishermitian::Bool, tags::ITensors.TagSets.GenericTagSet{…}, lefttags::Nothing, righttags::Nothing, plev::Nothing, leftplev::Nothing, rightplev::Nothing)
    @ ITensors ~/.julia/packages/ITensors/Xs85I/src/tensor_operations/matrix_decomposition.jl:383
  [4] contract(::NDTensors.BackendSelection.Algorithm{…}, A::MPO, ψ::MPS; cutoff::Float64, maxdim::Int64, mindim::Int64, normalize::Bool, kwargs::@Kwargs{})
    @ ITensorMPS ~/.julia/packages/ITensorMPS/AQY99/src/mpo.jl:743
  [5] contract(::NDTensors.BackendSelection.Algorithm{:densitymatrix, @NamedTuple{}}, A::MPO, ψ::MPS)
    @ ITensorMPS ~/.julia/packages/ITensorMPS/AQY99/src/mpo.jl:692
  [6] product(alg::NDTensors.BackendSelection.Algorithm{:densitymatrix, @NamedTuple{}}, A::MPO, ψ::MPS; kwargs::@Kwargs{})
    @ ITensorMPS ~/.julia/packages/ITensorMPS/AQY99/src/mpo.jl:605
  [7] product
    @ ~/.julia/packages/ITensorMPS/AQY99/src/mpo.jl:604 [inlined]
  [8] #apply#457
    @ ~/.julia/packages/ITensorMPS/AQY99/src/mpo.jl:601 [inlined]
  [9] product(A::MPO, ψ::MPS)
    @ ITensorMPS ~/.julia/packages/ITensorMPS/AQY99/src/mpo.jl:600
 [10] top-level scope
    @ REPL[7]:1

I don’t know what’s happening here… the MPO conserves the particle number so it should be applicable to the MPS. What am I doing wrong that triggers this error?

[miles]

I haven’t taken a very deep look at this yet, but I thought I’d try running your code and doing some checks of the inputs. I find that I can reproduce your same error. Something pops out to me though, which is that if I evaluate the QN flux of the xs MPS using @show flux(xs) I get that the flux is zero. However, when I measure expect(xs,"Ntot") I find that the expected value \langle N^{tot}_j \rangle is 2.0 for each site.

So I mention this to say something seems a bit off with the MPS used.

Where in the code do you actually initialize the MPS to be a certain physical state?

[corbett5]

Change it to

(t[coords.I..., 1, 1] != 0) && (newT[coords.I..., 1] = t[coords.I..., 1, 1])

Even if t[coords.I..., 1, 1] == 0 you’re still inserting a block with non-zero flux into the sparse tensor. I realized this after posting the code you linked but never got around to updating it, and now it’s too late. My bad.

[phx]

Thank you With your fix I don’t get that error anymore. The code I posted above now works.