# Issues with MPO generation with QN conservation **URL:** \ **Category:** ITensor Julia Questions\ **Tags:** mpo\ **Created:** [February 29, 2024, 5:09pm UTC](https://itensor.discourse.group/t/issues-with-mpo-generation-with-qn-conservation/1515 "2024-02-29T17:09:38Z")\ **Posts on this page:** 5\ **Page:** 1
**Author:** ![joshperry3110](https://avatars.discourse-cdn.com/v4/letter/j/839c29/32.png) [@joshperry3110](https://itensor.discourse.group/u/joshperry3110)\ **Post date:** [February 29, 2024, 5:09pm UTC](https://itensor.discourse.group/t/issues-with-mpo-generation-with-qn-conservation/1515/1 "2024-02-29T17:09:38Z")
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: ```Julia 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: ```Julia 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. ---
**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\ **Post date:** [February 29, 2024, 5:52pm UTC](https://itensor.discourse.group/t/issues-with-mpo-generation-with-qn-conservation/1515/2 "2024-02-29T17:52:57Z")
Thanks for asking. The issue you’re running into is that, at present, ITensor isn’t quite flexible enough to allow mixing QN-conserving tensors with regular, dense tensors. So the way you are using the code, it is trying to make dense tensors for the Boson sites (since for those you do not specify any QN’s should be conserved), while making block-sparse, QN-conserving tensors for the fermion sites. A solution would be to also say `conserve_qns=true` for the Boson sites. Can you try that and see if it works for you? Or is your Hamiltonian such that the total boson number is not conserved? Then I may be able to suggest a different approach. ---
**Author:** ![joshperry3110](https://avatars.discourse-cdn.com/v4/letter/j/839c29/32.png) [@joshperry3110](https://itensor.discourse.group/u/joshperry3110)\ **Post date:** [February 29, 2024, 7:58pm UTC](https://itensor.discourse.group/t/issues-with-mpo-generation-with-qn-conservation/1515/3 "2024-02-29T19:58:57Z")
Unfortunately the boson number is not conserved within my model. I would appreciate any insight into an alternative approach. Thanks ---
**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\ **Post date:** [February 29, 2024, 9:17pm UTC](https://itensor.discourse.group/t/issues-with-mpo-generation-with-qn-conservation/1515/4 "2024-02-29T21:17:44Z")
I see, so here is the trick you can use instead then. Each time you need to define one of your Boson sites, you can do it like this: ```julia sites[1] = Index([QN()=>N_max],"Boson,n=1") ``` or for other sites you can change `"n=1"` in the tag to `"n=7"` or `"n=$j"` for some site number `j`. The main thing is to replace the dimension so that it is no longer a number like `3` or `N_max` but instead now it is a vector with a single entry `[QN()=>N_max]`. In words, what that vector means is something like “there is one subspace with a dimension of N\_max whose quantum number is QN() namely the trivial quantum number”. So it’s like putting in “fake” quantum numbers so that the tensor will be internally a ‘quantum number conserving’ tensor but will actually just have a single block so will actually be equal to a dense tensor. Then you should be able to use it just like before except now it should be mixable with other actually sparse quantum number tensors like the ones associated to your fermion sites. ---
**Author:** ![joshperry3110](https://avatars.discourse-cdn.com/v4/letter/j/839c29/32.png) [@joshperry3110](https://itensor.discourse.group/u/joshperry3110)\ **Post date:** [March 1, 2024, 10:46am UTC](https://itensor.discourse.group/t/issues-with-mpo-generation-with-qn-conservation/1515/5 "2024-03-01T10:46:15Z")
Thank you for this, unfortunately I have now ran into another error as a result of this. I am trying to add together two MPO’s, one of which is the hermitian conjugate of the other. Previously I had this working but I believe the quantum number conserving sites has broken this for me. The error code is: ```julia LoadError: Attempting to contract IndexSet: ((dim=16|id=484|"Boson,n=1") 1: QN() => 16, (dim=2|id=987|"Link,l=1") 1: QN() => 1 2: QN() => 1, (dim=16|id=484|"Boson,n=1")' 1: QN() => 16) with IndexSet: ((dim=256|id=120|"Boson,n=1") 1: QN() => 256, (dim=16|id=484|"Boson,n=1")' 1: QN() => 16, (dim=16|id=484|"Boson,n=1") 1: QN() => 16) QN indices must have opposite direction to contract, but indices: (dim=16|id=484|"Boson,n=1") 1: QN() => 16 and: (dim=16|id=484|"Boson,n=1") 1: QN() => 16 do not have opposite directions. Stacktrace: [1] error(s::String) @ Base ./error.jl:35 [2] compute_contraction_labels(Ais::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, Bis::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}) @ ITensors ~/.julia/packages/ITensors/YVElE/src/indexset.jl:673 [3] _contract(A::NDTensors.BlockSparseTensor{ComplexF64, 3, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{ComplexF64, Vector{ComplexF64}, 3}}, B::NDTensors.Tensor{Number, 3, NDTensors.Combiner, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}}) @ ITensors ~/.julia/packages/ITensors/YVElE/src/tensor_operations/tensor_algebra.jl:2 [4] _contract(A::ITensor, B::ITensor) @ ITensors ~/.julia/packages/ITensors/YVElE/src/tensor_operations/tensor_algebra.jl:9 [5] contract(A::ITensor, B::ITensor) @ ITensors ~/.julia/packages/ITensors/YVElE/src/tensor_operations/tensor_algebra.jl:104 [6] * @ ~/.julia/packages/ITensors/YVElE/src/tensor_operations/tensor_algebra.jl:91 [inlined] [7] _broadcast_getindex_evalf @ ./broadcast.jl:683 [inlined] [8] _broadcast_getindex @ ./broadcast.jl:656 [inlined] [9] getindex @ ./broadcast.jl:610 [inlined] [10] macro expansion @ ./broadcast.jl:974 [inlined] [11] macro expansion @ ./simdloop.jl:77 [inlined] [12] copyto! @ ./broadcast.jl:973 [inlined] [13] copyto! @ ./broadcast.jl:926 [inlined] [14] copyto!(ψ::MPO, b::Base.Broadcast.Broadcasted{Base.Broadcast.Style{MPO}, Tuple{Base.OneTo{Int64}}, typeof(*), Tuple{MPO, Vector{ITensor}}}) @ ITensors ~/.julia/packages/ITensors/YVElE/src/mps/abstractmps.jl:2366 [15] materialize! @ ./broadcast.jl:884 [inlined] [16] materialize!(dest::MPO, bc::Base.Broadcast.Broadcasted{Base.Broadcast.Style{MPO}, Nothing, typeof(*), Tuple{MPO, Vector{ITensor}}}) @ Base.Broadcast ./broadcast.jl:881 [17] +(::NDTensors.Algorithm{:densitymatrix}, ::MPO, ::Vararg{MPO}; cutoff::Float64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @ ITensors ~/.julia/packages/ITensors/YVElE/src/mps/abstractmps.jl:210 [18] +(::NDTensors.Algorithm{:densitymatrix}, ::MPO, ::MPO) @ ITensors ~/.julia/packages/ITensors/YVElE/src/mps/abstractmps.jl:1380 [19] +(::MPO, ::Vararg{MPO}; alg::NDTensors.Algorithm{:densitymatrix}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @ ITensors ~/.julia/packages/ITensors/YVElE/src/mps/abstractmps.jl:1548 [20] + @ ~/.julia/packages/ITensors/YVElE/src/mps/abstractmps.jl:1547 [inlined] ``` This comes from the lines: ```julia X = dag(Kinetic_MPO) swapprime!(X, 0 => 1) H1 = X + Kinetic_MPO ``` Specifically the final line there. Is there a way to swap the direction of the index? Thank you so much