Hello,

I’m currently trying to define the sites for a boson-fermion model with some coupling. I am using one boson site and then a chain of L-1 fermions, and am trying to conserve the overall fermion number.

Previously before I added the QN conservation the code ran. However when I add conserve_qns=true to the fermion sites I get an error. Here is my code:

```
using ITensors
function Kinetic_MPO_Generator(L, sites, t_h)
op_kin = OpSum()
op_kin += "Id", 1
for j in 2:L-1
op_kin += "Cdag", j, "C", j+1
end
op_kin += "Cdag", L, "C", 2 #PBC
op_kin = -t_h * op_kin
Kinetic_MPO = MPO(op_kin, sites)
return Kinetic_MPO
end
L = 5
N_max = 8
t_h = 1
sites = Vector{Index}(undef, L)
sites[1] = siteind("Boson", dim=N_max)
for j=2:L
sites[j] = siteind("Fermion"; conserve_qns=true)
end
H = Kinetic_MPO_Generator(L, sites, t_h)
```

The line generating the error is the Kinetic_MPO = MPO(op_kin, sites) in the function. The error is as follows:

```
ERROR: LoadError: Can't contract tensor of storage type NDTensors.Dense{Float64, Vector{Float64}} with tensor of storage type NDTensors.BlockSparse{Float64, Vector{Float64}, 2}.
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] contract(::Type{SimpleTraits.Not{NDTensors.CanContract{NDTensors.DenseTensor{Float64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, NDTensors.BlockSparseTensor{Float64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 2}}}}}, tensor1::NDTensors.DenseTensor{Float64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, labels_tensor1::Tuple{Int64, Int64}, tensor2::NDTensors.BlockSparseTensor{Float64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 2}}, labels_tensor2::Tuple{Int64, Int64})
@ NDTensors ~/.julia/packages/NDTensors/kR4oQ/src/tensoralgebra/generic_tensor_operations.jl:82
[3] contract
@ ~/.julia/packages/SimpleTraits/l1ZsK/src/SimpleTraits.jl:331 [inlined]
[4] _contract(A::NDTensors.DenseTensor{Float64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, B::NDTensors.BlockSparseTensor{Float64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 2}})
@ ITensors ~/.julia/packages/ITensors/YVElE/src/tensor_operations/tensor_algebra.jl:3
[5] _contract(A::ITensor, B::ITensor)
@ ITensors ~/.julia/packages/ITensors/YVElE/src/tensor_operations/tensor_algebra.jl:9
[6] contract(A::ITensor, B::ITensor)
@ ITensors ~/.julia/packages/ITensors/YVElE/src/tensor_operations/tensor_algebra.jl:104
[7] *
@ ~/.julia/packages/ITensors/YVElE/src/tensor_operations/tensor_algebra.jl:91 [inlined]
[8] svdMPO(ValType::Type{Float64}, os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index}; mindim::Int64, maxdim::Int64, cutoff::Float64)
@ ITensors ~/.julia/packages/ITensors/YVElE/src/physics/autompo/opsum_to_mpo.jl:116
[9] svdMPO
@ ~/.julia/packages/ITensors/YVElE/src/physics/autompo/opsum_to_mpo.jl:3 [inlined]
[10] svdMPO(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ITensors ~/.julia/packages/ITensors/YVElE/src/physics/autompo/opsum_to_mpo.jl:145
[11] svdMPO
@ ~/.julia/packages/ITensors/YVElE/src/physics/autompo/opsum_to_mpo.jl:142 [inlined]
[12] MPO(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index}; splitblocks::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ITensors ~/.julia/packages/ITensors/YVElE/src/physics/autompo/opsum_to_mpo_generic.jl:302
[13] MPO
@ ~/.julia/packages/ITensors/YVElE/src/physics/autompo/opsum_to_mpo_generic.jl:292 [inlined]
[14] Kinetic_MPO_Generator(L::Int64, sites::Vector{Index}, t_h::Int64)
```

If I remove the conserve_qns=true part from the function, the code runs fine.

Any help with this issue would be greatly appreciated. Thank you.