Mixed-site (Fermion,boson) model: Error calling MPO(ampo,sites)

Hi everyone,

calling the MPO method when constructing Hamiltonian with both fermions and bosons throws the error

ERROR: MethodError: no method matching op(::OpName{:F}, ::SiteType{Qudit}, ::Int64)
Closest candidates are:
  op(::Function, ::Any...; kwargs...) at /Users/rudolfsmorka/.julia/packages/ITensors/5CAqA/src/physics/sitetype.jl:458
  op(::Vector{var"#s584"} where var"#s584"<:Index, ::Any, ::Integer...; kwargs...) at /Users/rudolfsmorka/.julia/packages/ITensors/5CAqA/src/physics/sitetype.jl:442
  op(::Vector{var"#s586"} where var"#s586"<:Index, ::Any, ::Integer, ::NamedTuple) at /Users/rudolfsmorka/.julia/packages/ITensors/5CAqA/src/physics/sitetype.jl:450
  ...
Stacktrace:
  [1] op(on::OpName{:F}, st::SiteType{Boson}, ds::Int64; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors ~/.julia/packages/ITensors/5CAqA/src/physics/site_types/boson.jl:24
  [2] op(on::OpName{:F}, st::SiteType{Boson}, ds::Int64)
    @ ITensors ~/.julia/packages/ITensors/5CAqA/src/physics/site_types/boson.jl:24
  [3] op(::OpName{:F}, ::SiteType{Boson}, ::Index{Int64}; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors ~/.julia/packages/ITensors/5CAqA/src/physics/site_types/boson.jl:32
  [4] op(::OpName{:F}, ::SiteType{Boson}, ::Index{Int64})
    @ ITensors ~/.julia/packages/ITensors/5CAqA/src/physics/site_types/boson.jl:30
  [5] op(name::String, s::Index{Int64}; adjoint::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors ~/.julia/packages/ITensors/5CAqA/src/physics/sitetype.jl:289
  [6] op
    @ ~/.julia/packages/ITensors/5CAqA/src/physics/sitetype.jl:226 [inlined]
  [7] #op#1018
    @ ~/.julia/packages/ITensors/5CAqA/src/physics/sitetype.jl:410 [inlined]
  [8] op(s::Index{Int64}, opname::String)
    @ ITensors ~/.julia/packages/ITensors/5CAqA/src/physics/sitetype.jl:410
  [9] computeSiteProd(sites::Vector{Index{Int64}}, ops::Prod{Op})
    @ ITensors ~/.julia/packages/ITensors/5CAqA/src/physics/autompo/opsum_to_mpo_generic.jl:80
 [10] svdMPO(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index{Int64}}; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors ~/.julia/packages/ITensors/5CAqA/src/physics/autompo/opsum_to_mpo.jl:119
 [11] svdMPO
    @ ~/.julia/packages/ITensors/5CAqA/src/physics/autompo/opsum_to_mpo.jl:2 [inlined]
 [12] MPO(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index{Int64}}; splitblocks::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors ~/.julia/packages/ITensors/5CAqA/src/physics/autompo/opsum_to_mpo_generic.jl:248
 [13] MPO
    @ ~/.julia/packages/ITensors/5CAqA/src/physics/autompo/opsum_to_mpo_generic.jl:239 [inlined]
 [14] HamiltonianMPO(sites::Vector{Index{Int64}}, N::Int64, Eon::Float64, t0::Float64, gPh::Float64)
    @ Main ./REPL[4]:15
 [15] top-level scope
    @ REPL[17]:1

of which I am not sure how to handle.
The cpp equivalent of this code does not throw this error. I checked that the site indices match the particle type, and they do. However, the error message suggests, that they do not match (if it is correct to interpret “OpName{:F}” as a Fermion-operator name), which I find strange (without prior knowledge on how the Hamiltonian MPO is constructed internally in iTensor).

The example I used is a standard Holstein-model

H =\sum_{i} \varepsilon_i^{\phantom{\dagger}} c_i^\dagger c_{i}^{\phantom{\dagger}} -t_0\sum_{i} (c_i^\dagger c_{i+1}^{\phantom{\dagger}} +\mathrm{h.c.})+ g\sum_i n_i (a_i^\dagger + a_i^{\phantom{\dagger}}).

The code I use is:

using ITensors
function HamiltonianMPO(sites,N,Eon,t0,gPh)
        ampo = OpSum()
        for i in 1:2:N
                ampo += (Eon,"Cdag",i,"C",i)
        end
        for i in 1:2:N-3
                ampo += (-t0,"Cdag",i,"C",i+2)
                ampo += (-t0,"Cdag",i+2,"C",i)
        end
        for i in 1:2:N-1
                ampo += (gPh,"n",i,"Adag",i+1)
                ampo += (gPh,"n",i,"A",i+1)
        end
        println(ampo)
        return MPO(ampo,sites)
end
fS = 3
N = 2*fS
dt = 0.01
tend = 5
Eon = 0.0
t0 = 1.0
gPh = 0.0

sites = siteinds(n->isodd(n) ? "Fermion" : "Boson",N)
H = HamiltonianMPO(sites,N,Eon,t0,gPh)

Thanks very much!

Hi Rudi,
Sorry you ran into this issue. It is a bug and I just submitted a pull request to fix it. So it will be fixed in the next version of ITensors.jl, but in the meantime, just add the following line of code to the top of your code and it ought to fix it for you right away:

ITensors.op(on::OpName"F", st::SiteType"Qudit", ds::Int...) = op(OpName"Id"(), st, ds...)

Please let me know if that doesn’t work of course.