four point correlation at same site

hello, I think I meet some troubles at four point MPO, this is a simple example code

using ITensors

let
    L = 4
    n_max = 2
    Np = 2
    sites = siteinds("Boson", L; dim=n_max, conserve_qns=true)
    state = ["0", "1", "0", "1"]
    psi = MPS(sites, state)
    sites = siteinds(psi)
    os = OpSum()
    os += "Adag", 1, "A", 2, "Adag", 1, "A", 2
    @show os
    op = MPO(os, sites)
    corr = inner(psi', op, psi)
    @show corr
end

for n_max > 2, it works well
for n_max = 2, os is not at same site, for example, os += “Adag”, 1, “A”, 2, “Adag”, 2, “A”, 3, it works well
but for n_max = 2, os is at same site, for example, 1, 2, 1, 2, it can not work and throw a error:
ERROR: LoadError: Could not find block of QNIndex with matching QN
how can I repair it? thanks

I think, for this kind of thing, it’s best to keep the site indices (such as 1,2,1,2) from matching each other and then work out the cases where they do match by hand on paper, then write them as a separate piece of code.

I’m fairly sure the reason for the error you are getting is that the operator is turning out to be zero (because a^\dagger_i a^\dagger_i = 0 for dimension 2 bosons), but instead of the error reporting it as zero, it is saying that there are no blocks (i.e. no non-zero blocks, meaning the only blocks are zero).

Thanks , by your reply I’m absolutely sure where my mistake is, and I solved my problem by convert the QN psi to dense(psi)., it works well.

@miles should we try to catch this kind of thing internally? I think it’s nice if it “just works” even if it is zero.

Good point – I’ll look into it and create an issue. Glancing at the code just now (which is inside the qn_svdMPO OpSum backend), it’s a somewhat tricky thing because it may require distinguishing inputs which give zero, as above, because they go outside the Hilbert space altogether, from inputs which do not respect the total QN flux of the operator, which should be treated as an erroneous input.