I am new to VUMPS, and I would like to learn how to use it as it can help me gain information about an infinite-sized chain following the Bose-Hubbard Hamiltonian.
I used the Extended Hubbard model example script available in the repository and modified it to work for a bosonic chain.
using ITensors, ITensorMPS
using ITensorInfiniteMPS
base_path = joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src")
src_files = ["vumps_subspace_expansion.jl", "entropy.jl"]
for f in src_files
include(joinpath(base_path, f))
end
#Bose-Hubbard model
function ITensorInfiniteMPS.unit_cell_terms(::Model"BH_chain"; J, U)
ampo = OpSum()
ampo += (-J,"Adag",1,"A",2)
ampo += (-J,"Adag",2,"A",1)
#on-site interactions
for j=1:2
ampo +=(0.5, U,"N",j,"N",j)
ampo +=(-0.5, U,"N",j)
end
return ampo
end
##############################################################################
# VUMPS parameters
#
maxdim = 50 # Maximum bond dimension
cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension
max_vumps_iters = 200 # Maximum number of iterations of the VUMPS algorithm at each bond dimension
vumps_tol = 1e-5
outer_iters = 5 # Number of times to increase the bond dimension
localham_type = MPO # or ITensor
conserve_qns = true
eager = true
model_params = (J=1.0, U=2.0)
##############################################################################
# CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
#
N = 2 # Unit cell size
@show N
@show localham_type
#initial guess: one boson per site (full filling)
initstate(n) = "1"
N_occ = 6
s = infsiteinds("Boson", N; initstate, conserve_qns=true, conserve_number=true, dim=N_occ+1)
ψ = InfMPS(s, initstate)
model = Model"BH_chain"()
@show model, model_params
# Form the Hamiltonian
H = InfiniteSum{localham_type}(model, s; model_params...)
# Check translational invariance
println("\nCheck translational invariance of initial infinite MPS")
@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))
outputlevel = 1
vumps_kwargs = (tol=vumps_tol, maxiter=max_vumps_iters, outputlevel, eager)
subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
# For now, to increase the bond dimension you must alternate
# between steps of VUMPS and subspace expansion (which outputs
# a new state that is equal to the original state but with
# a larger bond dimension)
println("\nRun VUMPS on initial product state, unit cell size $N")
ψ = vumps_subspace_expansion(H, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs)
# Check translational invariance
println("\nCheck translational invariance of optimized infinite MPS")
@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))
function expect_two_site(ψ::InfiniteCanonicalMPS, h::ITensor, n1n2)
n1, n2 = n1n2
ϕ = ψ.AL[n1] * ψ.AL[n2] * ψ.C[n2]
return (noprime(ϕ * h) * dag(ϕ))[]
end
function expect_two_site(ψ::InfiniteCanonicalMPS, h::MPO, n1n2)
return expect_two_site(ψ, prod(h), n1n2)
end
function expect_two_site(ψ::MPS, h::ITensor, n1n2)
n1, n2 = n1n2
ψ = orthogonalize(ψ, n1)
ϕ = ψ[n1] * ψ[n2]
return (noprime(ϕ * h) * dag(ϕ))[]
end
N_infinite = [expect(ψ, "N", n) for n in 1:N]
bs = [(1, 2), (2, 3)]
energy_infinite = map(b -> expect_two_site(ψ, H[b], b), bs)
#
# Compare to DMRG
#
Nfinite = 50
sfinite = siteinds("Boson", Nfinite; conserve_qns=true, conserve_number=true, dim=N_occ+1)
Hfinite = MPO(model, sfinite; model_params...)
ψfinite = random_mps(sfinite, initstate; linkdims=10)
println("\nQN sector of starting finite MPS")
@show flux(ψfinite)
nsweeps = 100
maxdims =
min.(maxdim, [2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 16, 16, 16, 16, 32, 32, 32, 32, 50])
@show maxdims
println("\nRun DMRG on $Nfinite sites")
energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite; nsweeps, maxdim=maxdims, cutoff)
println("\nEnergy density")
@show energy_finite_total / Nfinite
nfinite = Nfinite ÷ 2 - 1
bsfinite = [(nfinite, nfinite + 1), (nfinite + 1, nfinite + 2)]
hfinite(b) = ITensor(model, [sfinite[b[1]], sfinite[b[2]]]; model_params...)
energy_finite = map(b -> expect_two_site(ψfinite, hfinite(b), b), bsfinite)
N_finite = expect(ψfinite, "N")[nfinite:(nfinite + 1)]
#energy_exact = reference(model, Observable("energy"); model_params...)
corr_infinite = correlation_matrix(finite_mps(ψ, 1:10), "Adag", "A"; sites=2:11)
corr_finite = correlation_matrix(ψfinite, "Adag", "A"; sites=Int(Nfinite / 2):Int(Nfinite / 2 + 9))
S_finite = [entropy(ψfinite, b) for b in (Nfinite ÷ 2):(Nfinite ÷ 2 + N - 1)]
S_infinite = [entropy(ψ, b) for b in 1:N]
println("\nResults from VUMPS")
@show energy_infinite
@show N_infinite
@show corr_infinite
@show S_infinite
println("\nResults from DMRG")
@show energy_finite
@show N_finite
@show corr_finite
@show S_finite
nothing
If I were to run this code, the VUMPS and DMRG are giving same correlation matrices, population number and energy per cell, so I would trust that the routine is working.
However, I am a bit confused on my energies, as I get:
Results from VUMPS
energy_infinite = [-1.5124611797498109, 0.276393202250021]
Results from DMRG
energy_finite = [0.2763932022500227, -1.5124611797498133]
Energy density (from DMRG)
energy_finite_total / Nfinite = -0.5777126151554227
energy_infinite and energy_finite are similar, which is good, but they do not correspond to the resulting total energy / N . Is there something wrong? Or am I forgetting something important?