[WIP] Start simplifying the necessary OpSum definitions for Models

examples/development/finite_mps_to_infinite_mps.jl

@@ -0,0 +1,63 @@
+using ITensors
+using ITensorInfiniteMPS
+
+# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ
+function ising_opsum_infinite(N; J, h)
+  os = OpSum()
+  for n in 1:N
+    os -= J, "X", n, "X", n + 1
+  end
+  for n in 1:N
+    os -= h, "Z", n
+  end
+  return os
+end
+
+function ising_opsum_finite(N; J, h)
+  return ising_opsum_infinite(N - 1; J, h) + (-h, "Z", N)
+end
+
+function main(; N, J, h, nsites)
+  s = siteinds("S=1/2", N)
+
+  H = MPO(ising_opsum_finite(N; J=J, h=h), s)
+  ψ0 = randomMPS(s)
+
+  energy, ψ = dmrg(H, ψ0; nsweeps=10, cutoff=1e-10)
+  @show energy / N
+
+  n1 = N ÷ 2
+  n2 = n1 + 1
+
+  ψ = orthogonalize(ψ, n1)
+
+  println("\nFinite MPS energy on bond (n1, n2) = $((n1, n2))")
+
+  ϕ12 = ψ[n1] * ψ[n2]
+  ham12 = contract(MPO(ising_opsum_infinite(1; J=J, h=h), [s[n1], s[n2]]))
+  @show inner(ϕ12, apply(ham12, ϕ12))
+
+  # Get an infinite MPS that approximates
+  # the finite MPS in the bulk.
+  # `nsites` sets the number of sites in the unit
+  # cell of the infinite MPS.
+  ψ∞ = infinitemps_approx(ψ; nsites=nsites, nsweeps=10)
+
+  println("\nInfinite MPS energy on a few different bonds")
+  println("Infinite MPS has unit cell of size $nsites")
+
+  # Compute a few energies
+  for n1 in 1:3
+    n2 = n1 + 1
+    @show n1, n2
+    ϕ12 = ψ∞.AL[n1] * ψ∞.AL[n2] * ψ∞.C[n2]
+    s1 = siteind(ψ∞, n1)
+    s2 = siteind(ψ∞, n2)
+    ham12 = contract(MPO(ising_opsum_infinite(1; J=J, h=h), [s1, s2]))
+    @show inner(ϕ12, apply(ham12, ϕ12)) / norm(ϕ12)
+  end
+
+  return nothing
+end
+
+main(; N=100, J=1.0, h=1.5, nsites=1)

examples/development/orthogonalize_infinitemps.jl

@@ -0,0 +1,19 @@
+using ITensors
+using ITensorInfiniteMPS
+using Random
+
+Random.seed!(1234)
+
+# n-site unit cell
+nsites = 4
+s = siteinds("S=1/2", nsites; conserve_szparity=true)
+χ = 10
+@assert iseven(χ)
+space = (("SzParity", 1, 2) => χ ÷ 2) ⊕ (("SzParity", 0, 2) => χ ÷ 2)
+ψ = InfiniteMPS(ComplexF64, s; space=space)
+for n in 1:nsites
+  ψ[n] = randomITensor(inds(ψ[n]))
+end
+
+ψ = orthogonalize(ψ, :)
+@show norm(contract(ψ.AL[1:nsites]) * ψ.C[nsites] - ψ.C[0] * contract(ψ.AR[1:nsites]))

examples/vumps/transfer_matrix_spectrum.jl

@@ -94,7 +94,7 @@ nfinite = Nfinite ÷ 2
 orthogonalize!(ψfinite, nfinite)
 hnfinite = ITensor(model, sfinite[nfinite], sfinite[nfinite + 1]; model_params...)
 energy_finite = energy_local(ψfinite[nfinite], ψfinite[nfinite + 1], hnfinite)
-energy_infinite = energy_local(ψ.AL[1], ψ.AL[2] * ψ.C[2], H[(1, 2)])
+energy_infinite = energy_local(ψ.AL[1], ψ.AL[2] * ψ.C[2], contract(H[(1, 2)]))
 @show energy_finite, energy_infinite
 @show abs(energy_finite - energy_infinite)
 
@@ -115,7 +115,7 @@ Sz2_infinite = expect(ψ.AL[2] * ψ.C[2], "Sz")
 #
 
 using Arpack
-using KrylovKit
+using KrylovKit: eigsolve
 using LinearAlgebra
 
 T = TransferMatrix(ψ.AL)
@@ -138,7 +138,7 @@ tol = 1e-10
 neigs = length(v⃗ᴿ)
 
 # Normalize the vectors
-N⃗ = [(translatecell(v⃗ᴸ[n], 1) * v⃗ᴿ[n])[] for n in 1:neigs]
+N⃗ = [(translatecelltags(v⃗ᴸ[n], 1) * v⃗ᴿ[n])[] for n in 1:neigs]
 
 v⃗ᴿ ./= sqrt.(N⃗)
 v⃗ᴸ ./= sqrt.(N⃗)
@@ -152,18 +152,19 @@ function proj(v⃗ᴸ, v⃗ᴿ, n)
   Lⁿ = v⃗ᴸ[n]
   Rⁿ = v⃗ᴿ[n]
   return ITensorMap(
-    [translatecell(Lⁿ, 1), translatecell(Rⁿ, -1)]; input_inds=inds(Rⁿ), output_inds=inds(Lⁿ)
+    [translatecelltags(Lⁿ, 1), translatecelltags(Rⁿ, -1)]; input_inds=inds(Rⁿ), output_inds=inds(Lⁿ)
   )
 end
 
 P⃗ = [proj(v⃗ᴸ, v⃗ᴿ, n) for n in 1:neigs]
 T⁻P = T - sum(P⃗)
 
-#vⁱᴿ² = vⁱᴿ - (translatecell(v⃗ᴸ[1], 1) * vⁱᴿ)[] / norm(v⃗ᴿ[1]) * v⃗ᴿ[1]
+#vⁱᴿ² = vⁱᴿ - (translatecelltags(v⃗ᴸ[1], 1) * vⁱᴿ)[] / norm(v⃗ᴿ[1]) * v⃗ᴿ[1]
 #@show norm(dag(vⁱᴿ²) * v⃗ᴿ[1])
 
-λ⃗ᴾᴿ, v⃗ᴾᴿ, right_info = eigsolve(T⁻P, vⁱᴿ, neigs, :LM; tol=tol)
-@show λ⃗ᴾᴿ
+# XXX: Error with mismatched QN index arrows
+# λ⃗ᴾᴿ, v⃗ᴾᴿ, right_info = eigsolve(T⁻P, vⁱᴿ, neigs, :LM; tol=tol)
+# @show λ⃗ᴾᴿ
 
 vⁱᴿ⁻ᵈᵃᵗᵃ = vec(array(vⁱᴿ))
 λ⃗ᴿᴬ, v⃗ᴿ⁻ᵈᵃᵗᵃ = Arpack.eigs(T; v0=vⁱᴿ⁻ᵈᵃᵗᵃ, nev=neigs)

examples/vumps/vumps_heisenberg.jl

@@ -37,11 +37,11 @@ subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
 ψ = vumps(H, ψ; vumps_kwargs...)
 
 # Alternate steps of running VUMPS and increasing the bond dimension
-for _ in 1:outer_iters
+@time for _ in 1:outer_iters
   println("\nIncrease bond dimension")
-  global ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...)
+  global ψ = @time subspace_expansion(ψ, H; subspace_expansion_kwargs...)
   println("Run VUMPS with new bond dimension")
-  global ψ = vumps(H, ψ; vumps_kwargs...)
+  global ψ = @time vumps(H, ψ; vumps_kwargs...)
 end
 
 # Check translational invariance

examples/vumps/vumps_ising_extended.jl

@@ -7,15 +7,15 @@ using ITensorInfiniteMPS
 
 maxdim = 64 # Maximum bond dimension
 cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension
-max_vumps_iters = 100 # Maximum number of iterations of the VUMPS algorithm at each bond dimension
+max_vumps_iters = 10 # Maximum number of iterations of the VUMPS algorithm at each bond dimension
 vumps_tol = 1e-6
 outer_iters = 10 # Number of times to increase the bond dimension
 eager = true
 
 ##############################################################################
 # CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
 #
-N = 3# Number of sites in the unit cell
+N = 3 # Number of sites in the unit cell
 J = -1.0
 J₂ = -0.2
 h = 1.0;
@@ -53,17 +53,18 @@ sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
 Hfinite = MPO(model, sfinite; J=J, J₂=J₂, h=h)
 ψfinite = randomMPS(sfinite, initstate; linkdims=10)
 @show flux(ψfinite)
-sweeps = Sweeps(10)
-setmaxdim!(sweeps, maxdim)
-setcutoff!(sweeps, cutoff)
+dmrg_kwargs = (nsweeps=10, maxdim, cutoff)
+energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite; dmrg_kwargs...)
 energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps)
 
 nfinite = Nfinite ÷ 2
 Sz_finite = expect(ψfinite, "Sz")[nfinite:(nfinite + N - 1)]
 
+println("\nEnergy")
 @show energy_infinite
 @show energy_finite_total / Nfinite
 @show reference(model, Observable("energy"); J=J, h=h, J₂=J₂)
 
+println("\nSz")
 @show Sz
 @show Sz_finite

src/ITensorInfiniteMPS.jl

@@ -73,6 +73,7 @@ export Cell,
   reference,
   subspace_expansion,
   translatecell,
+  translatecelltags,
   translator,
   tdvp,
   vumps,

src/abstractinfinitemps.jl

@@ -206,13 +206,21 @@ end
 
 # TODO: return a Dictionary or IndexSetNetwork?
 function ITensors.siteinds(ψ::AbstractInfiniteMPS, r::AbstractRange)
-  return [siteinds(ψ, n) for n in r]
+  return [siteind(ψ, n) for n in r]
 end
 
-siterange(ψ::AbstractInfiniteMPS, c::Cell) = 1:(nsites(ψ) .+ (c.cell - 1))
+function siterange(ψ::AbstractInfiniteMPS, c::Cell)
+  #1:(nsites(ψ) .+ (c.cell - 1))
+  start = nsites(ψ) * (c.cell - 1) + 1
+  stop = start + nsites(ψ) - 1
+  return start:stop
+end
 
 ITensors.siteinds(ψ::AbstractInfiniteMPS, c::Cell) = siteinds(ψ, siterange(ψ, c))
 
+ITensors.siteinds(ψ::AbstractInfiniteMPS) = CelledVector(siteinds(ψ, Cell(1)), translator(ψ))
+infsiteinds(ψ::AbstractInfiniteMPS) = siteinds(ψ)
+
 Base.getindex(ψ::AbstractInfiniteMPS, r::UnitRange{Int}) = MPS([ψ[n] for n in r])
 
 struct UnitRangeToFunction{T<:Real,F<:Function}

src/celledvectors.jl

@@ -67,8 +67,9 @@ function translatecell(translator::Function, is::Union{<:Tuple,<:Vector}, n::Int
   return translatecell.(translator, is, n)
 end
 
-#Default behavior
-#translatecell(T::ITensor, n::Integer) = ITensors.setinds(T, translatecell(inds(T), n))
+# Default behavior
+translatecelltags(T::ITensor, n::Integer) = translatecell(translatecelltags, T, n)
+translatecelltags(T::ITensors.Indices, n::Integer) = translatecell(translatecelltags, T, n)
 #translatecell(T::MPO, n::Integer) = translatecell.(T, n)
 #translatecell(T::Matrix{ITensor}, n::Integer) = translatecell.(T, n)
 
@@ -99,11 +100,15 @@ function CelledVector{T}(v::Vector{T}, translator::Function) where {T}
 end
 
 """
-    celllength(cv::CelledVector)
+    cell_length(cv::CelledVector)
+
+    celllength(cv::CelledVector) # Deprecated
 
 The length of a unit cell of a CelledVector.
 """
-celllength(cv::CelledVector) = length(ITensors.data(cv))
+cell_length(cv::CelledVector) = length(ITensors.data(cv))
+
+celllength(cv::CelledVector) = cell_length(cv)
 
 # For compatibility with Base
 Base.size(cv::CelledVector) = size(ITensors.data(cv))