Hello ITensor community,

*Premise*

I would ultimately like to use the QN-conserving, block-sparse (possibly even not both at the same time) MPO representation of a long-range (LR) interacting Hamiltonian for the XXZ model;

$$

H = -J \sum_{n \neq m} \frac{1}{|n - m|^\alpha} \left[ \frac{1}{2}S^+*n S^-*m + \frac{1}{2}S^-*n S^+*m + \Delta S^z*n S^z*m \right]

$$

To efficiently represent such a Hamiltonian as an MPO, we use a sum of exponentially decaying interaction terms to mimic a power-law decay. This is a 'well known' technique as exponentially decaying interaction terms have literally zero extra cost in a (dense) MPO representation.

*Restricting my problem to minimal working example*

For the purpose of this question and for the minimal code example of the problem that I encounter, I will restrict myself to one single exponentially decaying interaction term in the Heisenberg S=1/2 model (here, the anisotropic coupling @@\Delta = 1@@). Although the bug is entirely independent of this choice.

*Quick background on exponentially decaying interaction range*

As, eg Schollwöck https://arxiv.org/abs/1008.3477 Eq. (188) explains, to construct an MPO with exponentially decaying interaction range, we simply have to introduce a Identity operator on the diagonal position corresponding to the interaction operator that is supposed to be coupled with LR interactions. In the bulk, an MPO has this matrix structure, when representing an interaction term like @@\sum_r J(r) S^z_i S^z_{i+r}@@:

$$

\hat{W}^{\mathrm{LRXXZ}} =\begin{bmatrix}

I & 0 & 0 \

S^z & \lambda I & 0 \

0 & J \lambda S^z & I

\end{bmatrix}

$$

**The actual problem:** Handmade QN conserving ITensor / MPS

So for the moment, I won't use a sum of exponentials which would trivially increase the MPO bond dimension by the amount of terms in the sum times the number of interaction terms in the Hamiltonian (3). For that reason, let us just use one weight for the single exponential (= 1.0) and one decay rate (@@\lambda = \exp(-1/\xi) =\exp(-1/4)@@).

I now wrote my little function that should craft myself a bulk MPO.

This function returns an MPO, with filled up ITensors. However, I seem to fail to find the mistake that I did since calling `dmrg(H, psi0, sweeps)`

fails.

Here would be a simple script that reproduces the Error:

```
using ITensors, Random
Random.seed!(421)
function build_long_range_HXXZ(sites, # ::Vector{Index{Vector{Pair{QN, Int64}}}}
Δ::T1,
α::T2,
onsite_couplings::Vector{<:Number};
kwargs...)::MPO where T1<:Number where T2<:Real
N = length(sites)
EType = eltype(union(Δ, α, onsite_couplings))
## ORDER OF INT_OPERATORS IS IMPORANT AND MUST MATCH THE QN()s of the linkdimensions
int_operators = iszero(Δ) ? [("S-", "S+"), ("S+", "S-") ] : [("Sz", "Sz"), ("S-", "S+"), ("S+", "S-")]
J = -1
int_couplings = iszero(Δ) ? J*[0.5, 0.5] : J*[0.5, 0.5, Δ]
onsite_operators = ["Sz"]
r_err::Float64 = get(kwargs, :r_err, 0.015)
r_abs::Float64 = get(kwargs, :a_err, 1e-3)
maxterms::Int64 = get(kwargs, :maxterms, 7)
# retrieve weights, bases and error
# weights, bases = get_weights_and_bases(N, Float64(α), r_err=r_err; maxterms=maxterms)
#### TEMPORARY annulation for demonstrative purposes
weights = [1.0]
bases = [exp(-1/4)]
#number of terms which are coupled via exponential decay
N_op = length(int_operators)
# check if h is actually zero everywhere
hasNoField = prod(map(iszero, onsite_couplings))
# if so, make N_op_onsite 0 manually and save the explicit 0 in the sparse matrix
N_op_onsite = hasNoField ? 0 : 1
# number of terms in sum of exponentials
M = length(weights)
if M > N÷2
@warn "Number of terms in sum of exponentials is exceeding half the system size!"
end
# link_dimension is given by produc of numbers of operators
# and numbers of terms in weights/bases (+2 for end and beginning)
d0 = 2
link_dimension = N_op*M + d0
sizeZeroBlock = N_op==2 ? d0 : d0 + M
# Magic part to be changed for different Hamitlonian#
linkindices = hasqns(sites) ? [Index([QN("Sz",0) => sizeZeroBlock, QN("Sz",-2) => M, QN("Sz",2) => M],"Link,l=$(n-1)") for n=1:N+1] : [Index(link_dimension,"Link,l=$(n-1)") for n=1:N+1] # QN index or "regular" index
# local_dim = dim(sites[2])
### n ∈[1,N] = 1
# W = Array{EType}(undef, link_dimension, link_dimension,local_dim,local_dim);
H = MPO(N)
for n=1:N
s = sites[n]
ll = dag(linkindices[n])
rl = linkindices[n+1]
# init empty ITensor with
H[n] = ITensor(EType, ll, rl, s', dag(s))
# add both Identities as netral elements in the MPS
H[n] += setelt(ll[1]) * setelt(rl[1]) * op(sites,"Id",n)
H[n] += setelt(ll[2]) * setelt(rl[2]) * op(sites,"Id",n)
# predefine indices m which count the M terms
ms = [0, 0, repeat(1:M,N_op)...]
# amd indices l which count N_op different interaction operators
ls = reduce(vcat, [0,0, [fill(i, M) for i in 1:N_op]...])
# loop over rest of columns 3:linkdim and fill out bulk MPS
for column in 3:dim(ll)
m = ms[column]
l = ls[column]
H[n] += weights[m] * bases[m] * int_couplings[l] * setelt(ll[2]) * setelt(rl[column]) * op(sites, int_operators[l][1],n) # α_m β_m Int_Op1, in the second row
H[n] += bases[m] * setelt(ll[column]) * setelt(rl[column]) * op(sites, "Id",n) # β_m Id, on the diagonal
H[n] += setelt(ll[column]) * setelt(rl[1]) * op(sites, int_operators[l][2],n) # Int_Op2, on the first column
# H[n][]
end
# on-site terms
if N_op_onsite > 0
for k = 1:N_op_onsite
H[n] += onsite_couplings[k] * setelt(ll[2]) * setelt(rl[1]) * op(sites, onsite_operators[k], n)
end
end
end
LE = setelt((linkindices[1][link_dimension]))
RE = setelt(dag(linkindices[N+1][1]))
H[1] *= LE # project out line vector because of OBC
H[N] *= RE # project out column vector because of OBC
return H
end
N = 12
Δ = 0.5
# α = 3.85
α = 10.0
h = 0.0
sites = siteinds("S=1/2", N; conserve_qns=true)
state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
psi0 = randomMPS(sites, state, 10)
sweeps = Sweeps(17)
mindim!(sweeps, 15, 15, 20, 40)
maxdim!(sweeps, 20, 20, 20, 30, 40, 60, 80, 100, 150, 400)
cutoff!(sweeps, 1e-10)
noise!(sweeps, 1e-5, 1e-5, 1e-5, 1e-6, 1e-6, 1e-7, 1e-7, 1e-7, 1e-8, 1e-8, 1e-8, 1e-8, 1e-8, 1e-9, 0.0)
H = build_long_range_HXXZ(sites, Δ, α, fill(h,N))
energy_gs, psi_gs = dmrg(H, psi0, sweeps; outputlevel=1)
```

This results in the error message:

```
julia> energy_gs, psi_gs = dmrg(H, psi0, sweeps; outputlevel=1)
ERROR: Eigen currently only supports block diagonal matrices.
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:33
[2] eigen(T::LinearAlgebra.Hermitian{Float64, ITensors.NDTensors.BlockSparseTensor{Float64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, ITensors.NDTensors.BlockSparse{Float64, Vector{Float64}, 2}}}; kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{10, Symbol}, NamedTuple{(:ishermitian, :which_decomp, :tags, :maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho, :normalize, :svd_alg), Tuple{Bool, Nothing, TagSet, Int64, Int64, Float64, ITensor, String, Bool, String}}})
@ ITensors.NDTensors ~/.julia/packages/ITensors/44bZA/src/NDTensors/blocksparse/linearalgebra.jl:203
[3] eigen(A::ITensor, Linds::Vector{Index{Vector{Pair{QN, Int64}}}}, Rinds::Vector{Index{Vector{Pair{QN, Int64}}}}; kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{10, Symbol}, NamedTuple{(:ishermitian, :which_decomp, :tags, :maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho, :normalize, :svd_alg), Tuple{Bool, Nothing, TagSet, Int64, Int64, Float64, ITensor, String, Bool, String}}})
@ ITensors ~/.julia/packages/ITensors/44bZA/src/decomp.jl:288
[4] factorize_eigen(A::ITensor, Linds::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}; kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{9, Symbol}, NamedTuple{(:which_decomp, :tags, :maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho, :normalize, :svd_alg), Tuple{Nothing, TagSet, Int64, Int64, Float64, ITensor, String, Bool, String}}})
@ ITensors ~/.julia/packages/ITensors/44bZA/src/decomp.jl:421
[5] factorize(A::ITensor, Linds::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}; kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{9, Symbol}, NamedTuple{(:which_decomp, :tags, :maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho, :normalize, :svd_alg), Tuple{Nothing, TagSet, Int64, Int64, Float64, ITensor, String, Bool, String}}})
@ ITensors ~/.julia/packages/ITensors/44bZA/src/decomp.jl:505
[6] replacebond!(M::MPS, b::Int64, phi::ITensor; kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{8, Symbol}, NamedTuple{(:maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho, :normalize, :which_decomp, :svd_alg), Tuple{Int64, Int64, Float64, ITensor, String, Bool, Nothing, String}}})
@ ITensors ~/.julia/packages/ITensors/44bZA/src/mps/mps.jl:453
[7] macro expansion
@ ~/.julia/packages/ITensors/44bZA/src/mps/dmrg.jl:250 [inlined]
[8] macro expansion
@ ~/.julia/packages/TimerOutputs/SSeq1/src/TimerOutput.jl:252 [inlined]
[9] macro expansion
@ ~/.julia/packages/ITensors/44bZA/src/mps/dmrg.jl:249 [inlined]
[10] macro expansion
@ ./timing.jl:287 [inlined]
[11] dmrg(PH::ProjMPO, psi0::MPS, sweeps::Sweeps; kwargs::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:outputlevel,), Tuple{Int64}}})
@ ITensors ~/.julia/packages/ITensors/44bZA/src/mps/dmrg.jl:188
[12] #dmrg#895
@ ~/.julia/packages/ITensors/44bZA/src/mps/dmrg.jl:47 [inlined]
[13] top-level scope
@ REPL[19]:1
```