Hi All,
I am trying to perform 2-D DMRG on a simple model (spinful free fermions) using a fractal “Hilbert curve” path as opposed to the more conventional “snake”. This strategy is proposed in [2507.11820v1] Fractal Path Strategies for Efficient 2D DMRG Simulations. The short description of the error is that combining S_z conservation with the choice of using a “Hilbert curve DMRG snake” leads to the LoadError: Eigen currently only supports block diagonal matrices. upon running the DMRG algorithm. If I use S_z conservation with a more conventional DMRG snake (the naive keyword in the code below) or if I use the Hilbert curve snake without S_z conservation, there is no error.
Given the above, I am sure that the error is not something silly like having OpSum terms that do not manifestly conserve the QN. I would have guessed that the long-range nature of the Hamiltonian (enforced by the Hilbert curve path) is responsible for this weird interaction with S_z conservation, but any 2D DMRG implementation should suffer from this issue.
Any advice is welcome, as well as any insight into how the QN functionality works and interacts with long-range interactions.
The code is provided at the very end in terms of functions; to run the computation, please use
nsweeps = 10, maxdim = 50;
conserve_sz = true #Toggle between false and true as needed
restart = false;
snaketype = "naive"; #Change to "hilbert" to test the fractal curve
nfl_spin_main(; Nx=4, Ny=4, t_SING, filling, nsweeps, maxdim, threaded_blocksparse=false, restart, conserve_sz, snaketype);
using ITensors, ITensorMPS, ITensorInfiniteMPS
using LinearAlgebra
using Random
using Strided
using HDF5
function FL_spinful(; Nx::Int, Ny::Int, t_SING=0.0, yperiodic::Bool=true,snaketype="naive")
N = Nx * Ny
lattice = square_lattice(Nx, Ny; yperiodic=yperiodic,snaketype)
opsum = OpSum()
for b in lattice
bondtype = bond_type(b,Ny)
tmat = [1 0; 0 1]
opsum -= t_SING*tmat[1,1], "Cdagup", b.s1, "Cup", b.s2
opsum -= t_SING*tmat[1,2], "Cdagup", b.s1, "Cdn", b.s2
opsum -= t_SING*tmat[2,1], "Cdagdn", b.s1, "Cup", b.s2
opsum -= t_SING*tmat[2,2], "Cdagdn", b.s1, "Cdn", b.s2
opsum -= conj(t_SING*tmat[1,1]), "Cdagup", b.s2, "Cup", b.s1
opsum -= conj(t_SING*tmat[2,1]), "Cdagup", b.s2, "Cdn", b.s1
opsum -= conj(t_SING*tmat[1,2]), "Cdagdn", b.s2, "Cup", b.s1
opsum -= conj(t_SING*tmat[2,2]), "Cdagdn", b.s2, "Cdn", b.s1
end
return opsum
end
function square_lattice_snaked(Nx::Int, Ny::Int; yperiodic=false, snaketype="naive")::Lattice
yperiodic = yperiodic && (Ny > 2)
N = Nx * Ny
Nbond = 2N - Ny + (yperiodic ? 0 : -Nx)
latt = Lattice(undef, Nbond)
b = 0
for n in 1:N
x, y = snake_to_xy(n, Nx, Ny, snaketype)
if x < Nx
latt[b += 1] = LatticeBond(n, xy_to_snake(x+1,y,Nx,Ny,snaketype), x, y, x + 1, y)
end
if Ny > 1
if y < Ny
latt[b += 1] = LatticeBond(n, xy_to_snake(x,y+1,Nx,Ny,snaketype), x, y, x, y + 1)
end
if yperiodic && y == 1
#latt[b += 1] = LatticeBond(n, n + Ny - 1, x, y, x, y + Ny - 1)
latt[b += 1] = LatticeBond(xy_to_snake(x,y+Ny-1,Nx,Ny,snaketype), n, x, y + Ny - 1, x, y)
end
end
end
return latt
end
And the Lattice and LatticeBond structs are as defined in ITensorMPS.jl/src/lattices/lattices.jl at e9650e5a2070690cd51ddc1f54644007189de537 · ITensor/ITensorMPS.jl · GitHub. The function to start the DMRG computation is below; the input snaketype is set to “naive” by default (how one usually does 2d DMRG) and one can set it to to “hilbert” when using the Hilbert space-filling curve.
function main(;
Nx::Int=8,
Ny::Int=4,
t_SING::Float64=1.0,
filling::Int=1,
filling_plus_x::Int=0,
maxdim::Int=300,
conserve_ky=false,
threaded_blocksparse=false,
nsweeps=15,
yperiodic=true,
random_init=true,
seed=1234,
restart = false,
conserve_sz = true,
writing=true,
n_preheating=0,
snaketype="naive"
)
# Helps make results reproducible when comparing
# sequential vs. threaded.
itensor_rng = Xoshiro()
Random.seed!(itensor_rng, seed)
@show Threads.nthreads()
# Disable other threading
BLAS.set_num_threads(1)
Strided.set_num_threads(1)
ITensors.enable_threaded_blocksparse(threaded_blocksparse)
@show ITensors.using_threaded_blocksparse()
N = Nx * Ny
cutoff = [1e-6]
noise = [1e-6, 1e-7, 1e-8, 0.0]
sites = siteinds("Electron", N; conserve_qns=true, conserve_sz)
os = NFL_spinful(; Nx, Ny, t,h,t_SING,U,U_NN, JH, yperiodic,snaketype)
# Create starting state with checkerboard
# pattern
state = map(CartesianIndices((Ny,Nx))) do I #to handle filling 1+x, I may need to move away from "do" syntax. See Ryan Levy's comment...
y=I[1]
x=I[2]
return (mod(x+y,filling)==0) ? (iseven(div(x+y,filling)) ? "Dn" : "Up") : "0"
end
display(state)
sz_append = conserve_sz ? "_szcons" : ""
restart_dirname = "nfl_spinful_NN_filling_"*string(filling)*"_t_"*string(t)*"_h_"*string(h)*"_U_"*string(U)*"_UNN_"*string(U_NN)*"_JH_"*string(JH)*"_t_SING_"*string(t_SING)*"_Nx_"*string(Nx)*"_Ny_"*string(Ny)*sz_append
restart_filename = restart_dirname*"/MPS_GS.h5"
psi0 = if random_init
random_mps(itensor_rng, sites, state; linkdims=2)
else
MPS(sites, state)
end
if restart && isdir(restart_dirname) && false
f = h5open(restart_filename,"r")
psi0 = read(f,"psi",MPS)
sites = siteinds(psi0)
psi0 = replace_siteinds(psi0, sites);
close(f)
end
if restart && isdir(restart_dirname)
psi0 = h5open(restart_filename,"r") do f
read(f,"psi",MPS)
end
sites = siteinds(psi0)
psi0 = replace_siteinds(psi0, sites);
#close(f)
#H = replace_siteinds(H, sites);
end
@show flux(psi0)
@show maxlinkdim(psi0)
H = MPO(os, sites)
# Number of structural nonzero elements in a bulk
# Hamiltonian MPO tensor
@show nnz(H[end ÷ 2])
@show nnzblocks(H[end ÷ 2])
@time energy, psi = @time dmrg(H, psi0; nsweeps, maxdim,mindim, cutoff, noise)
@show flux(psi)
@show maxlinkdim(psi)
There is an additional block of code for defining the functions xy_to_snake and snake_to_xy, included below:
# Stores generated Gilbert curves. First entry is a pair (Nx, Ny), and second is the map defining the curve
const _gilbertcache = IdDict{Tuple{Int,Int}, Vector{CartesianIndex{2}} }()
"""
gilbertindices(dims::Tuple{Int,Int}; majdim=dims[1] >= dims[2] ? 1 : 2)
Constructs a vector of `CartesianIndex` objects, orderd by a generalized Hilbert
space-filling, starting at `CartesianIndex(1,1)`. It will end at
`CartesianIndex(dims[1],1)` if `majdim==1`, or `CartesianIndex(1,dims[2])` if `majdim==2`
(or the closest feasible point).
"""
gilbertindices(dims::Tuple{Int,Int}; kwargs...) =
gilbertorder(CartesianIndices(dims); kwargs...)
"""
gilbertorder(mat::AbstractMatrix; majdim=size(mat,1) >= size(mat,2) ? 1 : 2)
Constructs a vector of the elements of `mat`, ordered by a generalized Hilbert
space-filling curve. The list will start at `mat[1,1]`, and end at `mat[end,1]` if
`majdim==1` or `mat[1,end]` if `majdim==2` (or the closest feasible point).
"""
function gilbertorder(mat::AbstractMatrix{T}; majdim=size(mat,1) >= size(mat,2) ? 1 : 2) where {T}
list = sizehint!(T[], length(mat))
if majdim == 1
append_gilbert!(list, mat)
else
append_gilbert!(list, permutedims(mat,(2,1)))
end
end
function append_gilbert!(list, mat::AbstractMatrix)
# 1 |* |
# | ) |
# a |v |
a,b = size(mat)
if a == 1 || b == 1
# single in one dimension
append!(list, mat)
elseif 2a > 3b
# long case: split into two
# +-----+
# 1 |* |
# || |
# a2|v |
# +-----+
# |* |
# || |
# a |v |
# +-----+
a2 = div(a,2)
if isodd(a2) && a > 2
a2 += 1
end
append_gilbert!(list, mat[1:a2,:])
append_gilbert!(list, mat[a2+1:a,:])
else
# standard case: split into three
# b2
# +---+---+
# 1 |*->|* |
# | || |
# a2| || |
# +---+| |
# | || |
# a |<-*|v |
# +---+----+
a2 = div(a,2)
b2 = div(b,2)
if isodd(b2) && b > 2
b2 += 1
end
append_gilbert!(list, permutedims(mat[1:a2,1:b2],(2,1)))
append_gilbert!(list, mat[:,b2+1:b])
append_gilbert!(list, permutedims(mat[a:-1:a2+1,b2:-1:1],(2,1)))
end
end
"""
linearindices(list::Vector{CartesianIndex{2}})
Construct an integer matrix `M` containing the integers `1:length(list)` such that
`M[list[i]] == i`.
"""
function linearindices(list::Vector{CartesianIndex{2}})
cmax = maximum(list)
L = zeros(Int,cmax.I)
for (i,c) in enumerate(list)
L[c] = i
end
return L
end
function xy_to_snake(x::Int,y::Int,Nx::Int,Ny::Int,snaketype::String)
#Maps 2d coordinates onto a 1d snake
#valid snake types: naive, naive_y, hilbert
#naive_y and hilbert MUST have periodic set to false
ind = 0
if snaketype == "naive"
ind = y+(x-1)*Ny
elseif snaketype == "naive_y"
ind = x+(y-1)*Ny
elseif snaketype == "hilbert"
#ind = xy_to_hilbert(x, y, Nx, Ny)
maplist = get!(_gilbertcache, (Nx, Ny)) do #Get or create the Hilbert map
list = gilbertindices((Nx,Ny))
end
#mapmatrix = GilbertCurves.linearindices(maplist) #fold it into a matrix for ease
mapmatrix = linearindices(maplist) #fold it into a matrix for ease
ind = mapmatrix[x,y]
else
throw("Not a valid snake type")
end
return ind
end
function snake_to_xy(n::Int,Nx::Int,Ny::Int,snaketype::String)
#Maps 1d snake onto a 2d coordinate
#valid snake types: naive, naive_y, hilbert
#naive_y and hilbert MUST have periodic set to false
x=0
y=0
if snaketype == "naive"
x = div(n - 1, Ny) + 1
y = mod(n - 1, Ny) + 1
elseif snaketype == "naive_y"
x = mod(n - 1, Nx) + 1
y = div(n - 1, Nx) + 1
elseif snaketype == "hilbert"
#coords = hilbert_to_xy(n, Nx, Ny)
#x = coords[1]
#y=coords[2]
maplist = get!(_gilbertcache, (Nx, Ny)) do #Get or create the Hilbert map
list = gilbertindices((Nx,Ny))
end
cartindex = maplist[n]
x=cartindex[1]
y=cartindex[2]
else
throw("Not a valid snake type")
end
return x,y
end