Issues with MPO generation with QN conservation

ITensor Julia Questions
1

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.

2

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.

1 Like
3

Unfortunately the boson number is not conserved within my model. I would appreciate any insight into an alternative approach.
Thanks

4

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:

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.

3 Likes
5

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:

LoadError: Attempting to contract IndexSet:

((dim=16|id=484|"Boson,n=1") <Out>
 1: QN() => 16, (dim=2|id=987|"Link,l=1") <Out>
 1: QN() => 1
 2: QN() => 1, (dim=16|id=484|"Boson,n=1")' <Out>
 1: QN() => 16)

with IndexSet:

((dim=256|id=120|"Boson,n=1") <In>
 1: QN() => 256, (dim=16|id=484|"Boson,n=1")' <Out>
 1: QN() => 16, (dim=16|id=484|"Boson,n=1") <Out>
 1: QN() => 16)

QN indices must have opposite direction to contract, but indices:

(dim=16|id=484|"Boson,n=1") <Out>
 1: QN() => 16

and:

(dim=16|id=484|"Boson,n=1") <Out>
 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:

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