Simplify the necessary OpSum definitions for Models (#63)

* New interface for models in terms of `ITensorInfiniteMPS.unit_cell_terms`, which should return an `OpSum` containing the terms of the Hamiltonian in the unit cell.

* Simplify the interface for `infsiteinds`, automatically shift the space of the site indices based on the initial state.

Project.toml

@@ -14,17 +14,19 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 
 [compat]
 Compat = "3, 4"
 HDF5 = "0.15, 0.16"
-ITensors = "0.3"
+ITensors = "0.3.18"
 Infinities = "0.1"
 IterTools = "1"
 KrylovKit = "0.5"
 OffsetArrays = "1"
 QuadGK = "2"
 Requires = "1"
+SplitApplyCombine = "1.2.2"
 julia = "1.6"
 
 [extras]

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/development/vumps/transfer_matrix_spectrum_arpack.jl

@@ -0,0 +1,193 @@
+using ITensors
+using ITensorInfiniteMPS
+
+##############################################################################
+# VUMPS parameters
+#
+
+maxdim = 10 # Maximum bond dimension
+cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension
+max_vumps_iters = 20 # Maximum number of iterations of the VUMPS algorithm at a fixed bond dimension
+tol = 1e-5 # Precision error tolerance for outer loop of VUMPS or TDVP
+outer_iters = 5 # Number of times to increase the bond dimension
+time_step = -Inf # -Inf corresponds to VUMPS
+solver_tol = (x -> x / 100) # Tolerance for the local solver (eigsolve in VUMPS and exponentiate in TDVP)
+multisite_update_alg = "parallel" # Choose between ["sequential", "parallel"]
+conserve_qns = false
+N = 2 # Number of sites in the unit cell (1 site unit cell is currently broken)
+
+# Parameters of the transverse field Ising model
+model_params = (J=1.0, h=0.9)
+
+##############################################################################
+# CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
+#
+
+model = Model("ising")
+
+initstate(n) = "↑"
+s = infsiteinds("S=1/2", N; initstate, conserve_szparity=conserve_qns)
+ψ = InfMPS(s, initstate)
+
+# Form the Hamiltonian
+H = InfiniteSum{MPO}(model, s; model_params...)
+
+# Check translational invariance
+@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))
+
+vumps_kwargs = (
+  tol=tol,
+  maxiter=max_vumps_iters,
+  solver_tol=solver_tol,
+  multisite_update_alg=multisite_update_alg,
+)
+subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
+#ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...)
+
+# Alternate steps of running VUMPS and increasing the bond dimension
+@time for _ in 1:outer_iters
+  println("\nIncrease bond dimension")
+  global ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...)
+  println("Run VUMPS with new bond dimension")
+  global ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...)
+end
+
+# Check translational invariance
+@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))
+
+#
+# Compare to DMRG
+#
+
+Nfinite = 100
+sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
+Hfinite = MPO(model, sfinite; model_params...)
+ψfinite = randomMPS(sfinite, initstate)
+@show flux(ψfinite)
+sweeps = Sweeps(10)
+setmaxdim!(sweeps, maxdim)
+setcutoff!(sweeps, cutoff)
+energy_finite_total, ψfinite = @time dmrg(Hfinite, ψfinite, sweeps)
+@show energy_finite_total / Nfinite
+
+function energy_local(ψ1, ψ2, h)
+  ϕ = ψ1 * ψ2
+  return (noprime(ϕ * h) * dag(ϕ))[]
+end
+
+function ITensors.expect(ψ, o)
+  return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[]
+end
+
+# Exact energy at criticality: 4/pi = 1.2732395447351628
+
+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], contract(H[(1, 2)]))
+@show energy_finite, energy_infinite
+@show abs(energy_finite - energy_infinite)
+
+energy_exact = reference(model, Observable("energy"); model_params...)
+@show energy_exact
+
+Sz1_finite = expect(ψfinite[nfinite], "Sz")
+orthogonalize!(ψfinite, nfinite + 1)
+Sz2_finite = expect(ψfinite[nfinite + 1], "Sz")
+Sz1_infinite = expect(ψ.AL[1] * ψ.C[1], "Sz")
+Sz2_infinite = expect(ψ.AL[2] * ψ.C[2], "Sz")
+
+@show Sz1_finite, Sz2_finite
+@show Sz1_infinite, Sz2_infinite
+
+##############################################################################
+# Compute eigenspace of the transfer matrix
+#
+
+using Arpack
+using KrylovKit: eigsolve
+using LinearAlgebra
+
+T = TransferMatrix(ψ.AL)
+Tᵀ = transpose(T)
+vⁱᴿ = randomITensor(dag(input_inds(T)))
+vⁱᴸ = randomITensor(dag(input_inds(Tᵀ)))
+
+neigs = 10
+tol = 1e-10
+λ⃗ᴿ, v⃗ᴿ, right_info = eigsolve(T, vⁱᴿ, neigs, :LM; tol=tol)
+λ⃗ᴸ, v⃗ᴸ, left_info = eigsolve(Tᵀ, vⁱᴸ, neigs, :LM; tol=tol)
+
+@show norm(T(v⃗ᴿ[1]) - λ⃗ᴿ[1] * v⃗ᴿ[1])
+@show norm(Tᵀ(v⃗ᴸ[1]) - λ⃗ᴸ[1] * v⃗ᴸ[1])
+
+@show λ⃗ᴿ
+@show λ⃗ᴸ
+@show flux.(v⃗ᴿ)
+
+neigs = min(length(v⃗ᴸ), length(v⃗ᴿ))
+v⃗ᴸ = v⃗ᴸ[1:neigs]
+v⃗ᴿ = v⃗ᴿ[1:neigs]
+
+# Normalize the vectors
+N⃗ = [(translatecelltags(v⃗ᴸ[n], 1) * v⃗ᴿ[n])[] for n in 1:neigs]
+
+v⃗ᴿ ./= sqrt.(N⃗)
+v⃗ᴸ ./= sqrt.(N⃗)
+
+# Form a second starting vector orthogonal to v⃗ᴿ[1]
+# This doesn't work. TODO: project out v⃗ᴿ[1], v⃗ᴸ[1] from T
+#λ⃗ᴿ², v⃗ᴿ², right_info_2 = eigsolve(T, vⁱᴿ², neigs, :LM; tol=tol)
+
+# Projector onto the n-th eigenstate
+function proj(v⃗ᴸ, v⃗ᴿ, n)
+  Lⁿ = v⃗ᴸ[n]
+  Rⁿ = v⃗ᴿ[n]
+  return ITensorMap(
+    [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ⁱᴿ - (translatecelltags(v⃗ᴸ[1], 1) * vⁱᴿ)[] / norm(v⃗ᴿ[1]) * v⃗ᴿ[1]
+#@show norm(dag(vⁱᴿ²) * v⃗ᴿ[1])
+
+# 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)
+
+## XXX: this is giving an error about trying to set the element of the wrong QN block for:
+## maxdim = 5
+## cutoff = 1e-12
+## max_vumps_iters = 10
+## outer_iters = 10
+## model_params = (J=1.0, h=0.8)
+##
+## v⃗ᴿᴬ = [itensor(v⃗ᴿ⁻ᵈᵃᵗᵃ[:, n], input_inds(T); tol=1e-4) for n in 1:length(λ⃗ᴿᴬ)]
+## @show flux.(v⃗ᴿᴬ)
+
+@show λ⃗ᴿᴬ
+
+# Full eigendecomposition
+
+Tfull = prod(T)
+DV = eigen(Tfull, input_inds(T), output_inds(T))
+
+@show norm(Tfull * DV.V - DV.Vt * DV.D)
+
+d = diag(array(DV.D))
+
+p = sortperm(d; by=abs, rev=true)
+@show p[1:neigs]
+@show d[p[1:neigs]]
+
+println("Error if ED with Arpack")
+@show d[p[1:neigs]] - λ⃗ᴿᴬ

examples/vumps/vumps_ising_noncontiguous.jl → examples/development/vumps/vumps_ising_noncontiguous.jl

@@ -18,13 +18,8 @@ N = 2# Number of sites in the unit cell
 J = 1.0
 h = 1.0;
 
-function space_shifted(q̃sz)
-  return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
-end
-
-space_ = fill(space_shifted(1), N);
-s = infsiteinds("S=1/2", N; space=space_)
 initstate(n) = "↑"
+s = infsiteinds("S=1/2", N; initstate)
 ψ = InfMPS(s, initstate);
 
 model = Model("ising");
@@ -70,7 +65,7 @@ Sz_finite = expect(ψfinite, "Sz")[nfinite:(nfinite + N - 1)]
 @show (
   energy_infinite,
   energy_finite_total / Nfinite,
-  reference(model, Observable("energy"); J=J, h=h, J₂=J₂),
+  reference(model, Observable("energy"); J=J, h=h),
 )
 
 @show (Sz, Sz_finite)

examples/vumps/vumps_localham.jl → examples/development/vumps/vumps_localham.jl

@@ -4,8 +4,6 @@ using Random
 
 Random.seed!(1234)
 
-include("models.jl")
-
 function _siteind(site_tag, n::Int; space)
   return addtags(Index(space, "Site,n=$n"), site_tag)
 end

examples/vumps/development/vumps_mpo.jl → examples/development/vumps/vumps_mpo.jl

@@ -10,35 +10,30 @@ maxiter = 20
 outer_iters = 3
 
 N = 2
-model = Model"ising"()
+model = Model("ising")
 
-function space_shifted(::Model"ising", q̃sz)
-  return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
-end
-
-space_ = fill(space_shifted(model, 1), N)
-s = infsiteinds("S=1/2", N; space=space_)
 initstate(n) = "↑"
+s = infsiteinds("S=1/2", N; initstate)
 ψ = InfMPS(s, initstate)
 
 model_params = (J=-1.0, h=0.9)
 vumps_kwargs = (multisite_update_alg="sequential", tol=tol, maxiter=maxiter)
 subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
 
-H = InfiniteITensorSum(model, s; model_params...)
+H = InfiniteSum{ITensor}(model, s; model_params...)
 # Alternate steps of running VUMPS and increasing the bond dimension
 ψ1 = vumps(H, ψ; vumps_kwargs...)
 for _ in 1:outer_iters
-  ψ1 = subspace_expansion(ψ1, H; subspace_expansion_kwargs...)
-  ψ1 = vumps(H, ψ1; vumps_kwargs...)
+  global ψ1 = subspace_expansion(ψ1, H; subspace_expansion_kwargs...)
+  global ψ1 = vumps(H, ψ1; vumps_kwargs...)
 end
 
 Hmpo = InfiniteMPOMatrix(model, s; model_params...)
 # Alternate steps of running VUMPS and increasing the bond dimension
 ψ2 = vumps(Hmpo, ψ; vumps_kwargs...)
 for _ in 1:outer_iters
-  ψ2 = subspace_expansion(ψ2, Hmpo; subspace_expansion_kwargs...)
-  ψ2 = vumps(Hmpo, ψ2; vumps_kwargs...)
+  global ψ2 = subspace_expansion(ψ2, Hmpo; subspace_expansion_kwargs...)
+  global ψ2 = vumps(Hmpo, ψ2; vumps_kwargs...)
 end
 
 SzSz = prod(op("Sz", s[1]) * op("Sz", s[2]))