Sites with three conserved QNs yield a wrong Hamiltonian

Hi,

I found a way to have two separately conserved QNs in the ITensor documentation. However, I got a problematic Hamiltonian when I did the same procedure for three conserved QNs. The same opsum yields two MPOs with different norms depending on whether I use three QNs.

I consider a bunch of decoupled Heisenberg chains. This model should have a U(1)^N symmetry. I can define the Hamiltonian in a standard way as follows:

using ITensorInfiniteMPS
function ITensorInfiniteMPS.unit_cell_terms(::Model"Heisenberg"; N::Int)
    opsum = OpSum()
    for i in 1:N
        opsum += 1/2, "S+", i, "S-", i+N
        opsum += 1/2, "S-", i, "S+", i+N
        opsum += "Sz", i, "Sz", i+N
    end
    return opsum
end

or in another way with extra zero terms:

function ITensorInfiniteMPS.unit_cell_terms(::Model"Hzeros"; N::Int, J4::Real = 0)
    opsum = OpSum()
    for i in 1:N
        opsum += 1/2, "S+", i, "S-", i+N
        opsum += 1/2, "S-", i, "S+", i+N
        opsum += "Sz", i, "Sz", i+N
    end
    for i in 1:N
		opsum += J4, "Sz", mod1(i-1,N)+N, "Sz", i, "Sz", i+N, "Sz", mod1(i+1,N)
    end
    return opsum
end

One can run the codes

using ITensorInfiniteMPS: opsum_finite
let L = 8, N = 3
    @show opsum_finite(Model("Heisenberg"), L*N; N=N)
    @show opsum_finite(Model("Hzeros"), L*N; N=N)
end

to see these two definitions only differ by terms with zero coefficient when the system is finite.

Nevertheless, the following codes

let L = 8, N = 3, conserve_qns = true
    custom_site(j) = addtags(siteind("S=1/2"; conserve_qns, qnname_sz = "Sz$(mod1(j,N))"), "n=$(j)")
    sitesN = map(custom_site, 1:N*L)
    sitesSz = siteinds("S=1/2", N*L; conserve_qns)

    h = MPO(opsum_finite(Model("Hzeros"), L*N; N=N), sitesN)
    H = MPO(opsum_finite(Model("Heisenberg"), L*N; N=N), sitesN);
    println("Norm difference in three-conserved QNs: ", norm(h - H))

    h = MPO(opsum_finite(Model("Hzeros"), L*N; N=N), sitesSz)
    H = MPO(opsum_finite(Model("Heisenberg"), L*N; N=N), sitesSz);
    println("Norm difference in a single QN: ", norm(h - H))
end

produce different norm differences from the two definitions. I also performed dmrg on the two Hamiltonians. When I use the definition with zeros in three-conserved QNs, the energy density per site does not agree with the analytical result e = 1/4 - \ln{2}.

Does anyone know why this discrepancy happens?