[Julia] DMRG with handmade H MPO

Question

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^zn S^zm \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

Comments

A little background on the MPO construction:

The bulk operator is a permuted version of the one in e.g. Schollwöck Eq. (184) review. This is because the operators of the XXZ model have a well defined flux and we wish to group all operators with the same flux of the conserved quantum number Sz:
- @@S^+_n S^-_m@@ has a net flux of 2 in ITensor units
- @@S^+_n S^-_m@@ has a net flux of -2 in ITensor units
- @@S^z_n S^z_m@@ has no flux (0) in ITensor (any) units

Consequently, we permute the last row vector in  Schollwöck Eq. (184) with the second row, and we permute the last and second to last column vectors with the second and third column. This way we grouped all locally Sz conserving terms in the uppermost block, all locally Sz changing terms (+/- 2) in the rest. Together with the insight of the exponentially decaying coupling, we have this bulk operator:

$$
\hat{W}^{\mathrm{LRXXZ}} =\begin{bmatrix}
I &  &  &  & \\
-hS^z & {I} & \lambda J\Delta S^z & \lambda (J/2)S^- & \lambda (J/2)S^+ \\
S^z & & \lambda I & & \\
S^+ & & & \lambda I & \\
S^- & & & & \lambda I
\end{bmatrix}
$$

---

EDIT: deleted "comment" for "answer"

Answer 1

So it turns out that making my problem presentable and revising a lot the code that I want to share made me find one mistake already!
Thanks for existing, ITensor community, otherwise I wouldn't have found it :D

The mistake snuck in when I copied the ITensors.jl template of manufactoring the bulk MPO. The mistake happens when I project out the first and left side of the bulk tensor to form open boundary conditions: Here I kept the last row vector of the bulk MPO. Correct would have been the second row vector, w.r.t. the bulk MPO presented in my first comment!

I fixed this and now it seems to work just fine!

Comments

Glad you already made progress so soon! Do you still have a question about this, or should we wait for you to post a separate one if so?

Miles