I’m writing a difficult to read code, just because I know the analytical result which is -0.274519. This value is meant to be compared with what is printed in the 4th column of this (a laptop is more than enough).
When i change Lanczos tolerance (line 45) i expect something to change. Default for Float64 is 1e-12. I even used 1e-21 and the result still changes on the sixth digit and I’m using Float64.
It seems that the best result (It’s an iterative variational algorithm, so it means it’s converging up to the analytical value) is reached with tol=1e-15. If i put less or more it’s more distant to the value I want it to reach.
Do not mistake what I wrote as tol_zero, that is something that let me supress numerical noise in the matrix and it saturates at 1e-14, which is what i used there.
You can change lanczos tolerance directly from line 45 and see how the result changes.
Lanczos is used on a sparse matrix and i use both eigenalues and eigenvectors.
using Printf
using SparseArrays
using Random
using LinearAlgebra
using BenchmarkTools
using KrylovKit
function IDMRG(tol_zeros::Real, Delta::Real, max_states::Integer, NIter::Integer, Energy::Vector, Energy_2::Vector, Energy_bond::Vector, lattice_size::Vector, Truncation_Error::Vector)
Random.seed!(1)
Sz = sparse(Float64[[0.5, 0] [0, -0.5]])
Sp = sparse(Float64[[0, 0] [1.0, 0]])
Sm = sparse(Float64[[0, 1.0] [0, 0]])
Ham = sparse(Float64[0 0; 0 0])
Block_Sp = Sp
Block_Sm = Sm
Block_Sz = Sz
Block_Ham = Ham
m::Int64 = 2 #hilbert size
temp_matr = Array{Float64}(undef, max_states*2, max_states)
for i = 1:NIter
lattice_size[i] = 2*i +2 #spatial row_size
#Adding 1 site
H_new = kron(Block_Ham, sparse(I,2,2)) .- Delta*kron(Block_Sz, Sz) .+ 0.5*(kron(Block_Sp, Sm) .+ kron(Block_Sm, Sp))
Block_Sz = kron(sparse(I,m,m), Sz) #definition of new Block_Sz
Block_Sp = kron(sparse(I,m,m), Sp) #definition of new Block_Sp
Block_Sm = kron(sparse(I,m,m), Sm) #definition of new Block_Sm
#SuperBlock
H_super = kron(H_new, sparse(I,m*2,m*2)) .+ kron(sparse(I,m*2,m*2), H_new) .- Delta*kron(Block_Sz, Block_Sz) .+ 0.5*(kron(Block_Sp, Block_Sm) .+ kron(Block_Sm, Block_Sp))
droptol!(H_super, tol_zeros)
#trunc index
m_next::Int64 = m*m*4 # hilbert space row_size SuperBlock
rho_dim::Int64 = m*2 #sqrt(m_next)
trunc_index::Int64 = min(rho_dim,max_states) #this is applied to density matrix diag
#Diagonalizing SuperBlock Lanczos
Energy_ground, evec_computed, info = eigsolve(H_super, rand(Float64, m_next), 1, :SR, Lanczos(krylovdim=10, tol=1e-15)) # maxiter = 100
if info.converged < 1
println("Lanczos did not converge in iteration: ", i)
exit()
end
Energy[i] = Energy_ground[1] #output of eigensolve is a vector of eigenvals
psi_ground = evec_computed[1]
Energy_2[i] = Energy[i]/lattice_size[i]
if i> 1
Energy_bond[i] = (Energy[i] - Energy[i-1]) / 2
else
Energy_bond[i] = Energy[i]/2
end
#Reduced Density Matrix
psi_matr = reshape(psi_ground, (rho_dim, rho_dim)) #psi matr is a view of psi allocated data
rho_matr = psi_matr * psi_matr' #a way of tracing rho_matrix on the right part
rho_matr_sym = Symmetric(rho_matr) #it seems this is indeed symmetric
#Diagonalizing Density Matrix, Trunc matrix
evals, evecs = eigen!(rho_matr_sym, sortby=-) #if I use eigen! i need to not use rho_matrx again
trunc_evals = @view evals[1:trunc_index]
trunc_matrx = @view evecs[:,1:trunc_index]
trunc_matrx .= ifelse.(abs.(trunc_matrx) .< tol_zeros, 0, trunc_matrx)
Truncation_Error[i] = 1 - sum(trunc_evals)
# Preallocation for mul!
temp = @view temp_matr[1:rho_dim,1:trunc_index]
#Higher dimension kron
mul!(temp, H_new, trunc_matrx)
Block_Ham = trunc_matrx' * temp
Block_Ham .= ifelse.(abs.(Block_Ham) .< tol_zeros, 0, Block_Ham)
mul!(temp, Block_Sz, trunc_matrx)
Block_Sz = trunc_matrx' * temp
Block_Sz .= ifelse.(abs.(Block_Sz) .< tol_zeros, 0, Block_Sz)
mul!(temp, Block_Sp, trunc_matrx)
Block_Sp = trunc_matrx' * temp
Block_Sp .= ifelse.(abs.(Block_Sp) .< tol_zeros, 0, Block_Sp)
mul!(temp, Block_Sm, trunc_matrx)
Block_Sm = trunc_matrx' * temp
Block_Sm .= ifelse.(abs.(Block_Sm) .< tol_zeros, 0, Block_Sm)
m = trunc_index # hilbert space row_size before next kron products
end
end
#Definition of starting parameters
tol_zeros = 1e-14
Delta::Float64 = 0.5
max_states::Int16 = 5
NIter::Int32 = 30 #Final lattice size: 2*NIter + 2
#output vectors
lattice_size = Vector{Int64}(undef, NIter)
Energy = Vector{Float64}(undef, NIter)
Energy_2 = Vector{Float64}(undef, NIter)
Energy_bond = Vector{Float64}(undef, NIter)
Truncation_Error = Vector{Float64}(undef, NIter)
IDMRG(tol_zeros, Delta, max_states, NIter, Energy, Energy_2, Energy_bond, lattice_size, Truncation_Error)
for i = 1:NIter
print(lpad(lattice_size[i], 2, "0"), " ")
@printf("| %.15e | %.15e | %.15e | %.15e | \n", Energy[i], Energy_bond[i], Energy_2[i], Truncation_Error[i])
end