# [Julia] DMRG with handmade H MPO

+1 vote
edited

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

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)
### n ∈[1,N] = 1
H = MPO(N)
for n=1:N
s = sites[n]

# 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) * setelt(rl) * op(sites,"Id",n)
H[n] += setelt(ll) * setelt(rl) * 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) * setelt(rl[column]) * op(sites, int_operators[l],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) * op(sites, int_operators[l],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) * setelt(rl) * op(sites, onsite_operators[k], n)
end
end

end
H *= 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:
 error(s::String)
@ Base ./error.jl:33
 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
 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
 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
 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
 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
 macro expansion
@ ~/.julia/packages/ITensors/44bZA/src/mps/dmrg.jl:250 [inlined]
 macro expansion
@ ~/.julia/packages/TimerOutputs/SSeq1/src/TimerOutput.jl:252 [inlined]
 macro expansion
@ ~/.julia/packages/ITensors/44bZA/src/mps/dmrg.jl:249 [inlined]
 macro expansion
@ ./timing.jl:287 [inlined]
 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
 #dmrg#895
@ ~/.julia/packages/ITensors/44bZA/src/mps/dmrg.jl:47 [inlined]
 top-level scope
@ REPL:1

commented by (200 points)
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}$$
commented by (200 points)
edited

+1 vote