[WIP] Start simplifying the necessary OpSum definitions for Models

examples/experimental/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/experimental/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/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/celledvectors.jl

@@ -99,11 +99,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))

src/infinitemps_approx.jl

@@ -1,4 +1,3 @@
-
 function siteind(ψ::InfiniteMPS, n::Integer)
   return uniqueind(ψ[n], ψ[n - 1], ψ[n + 1])
 end
@@ -37,7 +36,11 @@ function infinitemps_approx(
   A = InfiniteMPS(
     [settags(s[nrange[n]], addtags(site_tags, "n=$n")) for n in 1:length(nrange)]; space=χ
   )
-  randn!.(A)
+  # XXX: Why doesn't this work?
+  # randn!.(A)
+  for n in 1:N
+    A[n] = randomITensor(inds(A[n]))
+  end
   ψ∞ = orthogonalize(A, :)
 
   # Site indices of the infinite MPS

src/models/ising_extended.jl

@@ -1,48 +1,44 @@
-nrange(::Model"ising_extended") = 3
+# A vector of the terms associated with each site of the
+# unit cell.
 # H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ - J₂ Σⱼ XⱼZⱼ₊₁Xⱼ₊₂
-function ITensors.MPO(::Model"ising_extended", s; J=1.0, h=1.0, J₂=0.0)
-  N = length(s)
-  a = OpSum()
-  if J != 0
-    for n in 1:(N - 1)
-      a += -J, "X", n, "X", n + 1
-    end
-  end
-  if J₂ != 0
-    for n in 1:(N - 2)
-      a += -J₂, "X", n, "Z", n + 1, "X", n + 2
-    end
-  end
-  if h != 0
-    for n in 1:N
-      a += -h, "Z", n
-    end
-  end
-  return splitblocks(linkinds, MPO(a, s))
-end
-
-# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ - J₂ Σⱼ XⱼZⱼ₊₁Xⱼ₊₂
-function ITensors.OpSum(::Model"ising_extended", n1, n2; J=1.0, h=1.0, J₂=0.0)
+function unit_cell_terms(::Model"ising_extended"; J=1.0, h=1.0, J₂=0.0)
   opsum = OpSum()
-  if J != 0
-    opsum += -J / 2, "X", n1, "X", n2
-    opsum += -J / 2, "X", n2, "X", n2 + 1
-  end
-  if J₂ != 0
-    opsum += -J₂, "X", n1, "Z", n2, "X", n2 + 1
-  end
-  opsum += -h / 3, "Z", n1
-  opsum += -h / 3, "Z", n2
-  opsum += -h / 3, "Z", n2 + 1
-  return opsum
-end
-
-function ITensors.ITensor(::Model"ising_extended", s1, s2, s3; J=1.0, h=1.0, J₂=0.0)
-  opsum = OpSum(Model"ising_extended"(), 1, 2; J=J, h=h, J₂=J₂)
-  return prod(MPO(opsum, [s1, s2, s3]))
+  opsum -= J, "X", 1, "X", 2
+  opsum -= J₂, "X", 1, "Z", 2, "X", 3
+  opsum -= h, "Z", 1
+  return [opsum]
 end
 
 function reference(::Model"ising_extended", ::Observable"energy"; J=1.0, h=1.0, J₂=0.0)
   f(k) = sqrt((J * cos(k) + J₂ * cos(2k) - h)^2 + (J * sin(k) + J₂ * sin(2k))^2)
   return -1 / 2π * ITensorInfiniteMPS.∫(k -> f(k), -π, π)
 end
+
+#
+# Deprecated
+#
+
+# nrange(::Model"ising_extended") = 3
+
+# function ITensors.ITensor(::Model"ising_extended", s1::Index, s2::Index, s3::Index; J=1.0, h=1.0, J₂=0.0)
+#   opsum = OpSum(Model"ising_extended"(), 1, 2; J=J, h=h, J₂=J₂)
+#   return prod(MPO(opsum, [s1, s2, s3]))
+# end
+
+# # For finite Hamiltonian with open boundary conditions
+# # H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ - J₂ Σⱼ XⱼZⱼ₊₁Xⱼ₊₂
+# function opsum_finite(::Model"ising_extended", n; J=1.0, h=1.0, J₂=0.0)
+#   # TODO: define using opsum_infinite, filtering sites
+#   # that extend over the boundary.
+#   opsum = OpSum()
+#   for j in 1:n
+#     if j ≤ (n - 1)
+#       opsum += -J, "X", j, "X", j + 1
+#     end
+#     if j ≤ (n - 2)
+#       opsum += -J₂, "X", j, "Z", j + 1, "X", j + 2
+#     end
+#     opsum -= h, "Z", j
+#   end
+#   return opsum
+# end