Fix VUMPS for multi-site unit cells and nonlocal interactions (#31)

Author: LHerviou

examples/vumps/vumps_ising_extended.jl

@@ -0,0 +1,102 @@
+using ITensors
+using ITensorInfiniteMPS
+
+##############################################################################
+# VUMPS parameters
+#
+
+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
+vumps_tol = 1e-6
+outer_iters = 4 # Number of times to increase the bond dimension
+
+##############################################################################
+# CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
+#
+N = 3# Number of sites in the unit cell
+J = -1.0
+J₂ = -0.2
+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) = "↑"
+ψ = InfMPS(s, initstate);
+
+model = Model("ising_extended");
+H = InfiniteITensorSum(model, s; J=J, J₂=J₂, h=h);
+
+@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...));
+
+vumps_kwargs = (tol=vumps_tol, maxiter=max_vumps_iters)
+subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
+ψ_0 = vumps(H, ψ; vumps_kwargs...)
+
+for j in 1:outer_iters
+  println("\nIncrease bond dimension")
+  ψ_1 = subspace_expansion(ψ_0, H; subspace_expansion_kwargs...)
+  println("Run VUMPS with new bond dimension")
+  ψ_0 = vumps(H, ψ_1; vumps_kwargs...)
+end
+
+function ITensors.expect(ψ::InfiniteCanonicalMPS, o, n)
+  return (noprime(ψ.AL[n] * ψ.C[n] * op(o, s[n])) * dag(ψ.AL[n] * ψ.C[n]))[]
+end
+
+function ITensors.expect(ψ::InfiniteCanonicalMPS, h::ITensor)
+  l = linkinds(only, ψ.AL)
+  r = linkinds(only, ψ.AR)
+  s = siteinds(only, ψ)
+  δˢ(n) = δ(dag(s[n]), prime(s[n]))
+  δˡ(n) = δ(l[n], prime(dag(l[n])))
+  δʳ(n) = δ(dag(r[n]), prime(r[n]))
+  ψ′ = prime(dag(ψ))
+
+  ns = sort(findsites(ψ, h))
+  nrange = ns[end] - ns[1] + 1
+  idx = 2
+  temp_O = δˡ(ns[1] - 1) * ψ.AL[ns[1]] * h * ψ′.AL[ns[1]]
+  for n in (ns[1] + 1):(ns[1] + nrange - 1)
+    if n == ns[idx]
+      temp_O = temp_O * ψ.AL[n] * ψ′.AL[n]
+      idx += 1
+    else
+      temp_O = temp_O * ψ.AL[n] * δˢ(n) * ψ′.AL[n]
+    end
+  end
+  temp_O = temp_O * ψ.C[ns[end]] * δʳ(ns[end]) * ψ′.C[ns[end]]
+  return temp_O[]
+end
+
+function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteITensorSum)
+  return [expect(ψ, h[(j, j + 1)]) for j in 1:nsites(ψ)]
+end
+
+Sz = [expect(ψ_0, "Sz", n) for n in 1:N]
+energy_infinite = expect(ψ_0, H)
+
+Nfinite = 100
+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(15)
+setmaxdim!(sweeps, maxdim)
+setcutoff!(sweeps, cutoff)
+energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps)
+
+nfinite = Nfinite ÷ 2
+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₂),
+)
+
+@show (Sz, Sz_finite)

examples/vumps/vumps_ising_noncontiguous.jl

@@ -0,0 +1,105 @@
+using ITensors
+using ITensorInfiniteMPS
+
+##############################################################################
+# VUMPS parameters
+#
+
+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
+vumps_tol = 1e-6
+outer_iters = 4 # Number of times to increase the bond dimension
+
+##############################################################################
+# CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
+#
+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) = "↑"
+ψ = InfMPS(s, initstate);
+
+model = Model("ising");
+H = InfiniteITensorSum(model, s; J=J, h=h);
+#to test the case where the range is larger than the unit cell size
+for x in 1:N
+  H.data[x] = H.data[x] * delta(prime(s[x + 3]), dag(s[x + 3]))
+end
+
+@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...));
+
+vumps_kwargs = (tol=vumps_tol, maxiter=max_vumps_iters)
+subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
+ψ_0 = vumps(H, ψ; vumps_kwargs...)
+
+for j in 1:outer_iters
+  println("\nIncrease bond dimension")
+  ψ_1 = subspace_expansion(ψ_0, H; subspace_expansion_kwargs...)
+  println("Run VUMPS with new bond dimension")
+  ψ_0 = vumps(H, ψ_1; vumps_kwargs...)
+end
+
+function ITensors.expect(ψ::InfiniteCanonicalMPS, o, n)
+  return (noprime(ψ.AL[n] * ψ.C[n] * op(o, s[n])) * dag(ψ.AL[n] * ψ.C[n]))[]
+end
+
+function ITensors.expect(ψ::InfiniteCanonicalMPS, h::ITensor)
+  l = linkinds(only, ψ.AL)
+  r = linkinds(only, ψ.AR)
+  s = siteinds(only, ψ)
+  δˢ(n) = δ(dag(s[n]), prime(s[n]))
+  δˡ(n) = δ(l[n], prime(dag(l[n])))
+  δʳ(n) = δ(dag(r[n]), prime(r[n]))
+  ψ′ = prime(dag(ψ))
+
+  ns = sort(findsites(ψ, h))
+  nrange = ns[end] - ns[1] + 1
+  idx = 2
+  temp_O = δˡ(ns[1] - 1) * ψ.AL[ns[1]] * h * ψ′.AL[ns[1]]
+  for n in (ns[1] + 1):(ns[1] + nrange - 1)
+    if n == ns[idx]
+      temp_O = temp_O * ψ.AL[n] * ψ′.AL[n]
+      idx += 1
+    else
+      temp_O = temp_O * ψ.AL[n] * δˢ(n) * ψ′.AL[n]
+    end
+  end
+  temp_O = temp_O * ψ.C[ns[end]] * δʳ(ns[end]) * ψ′.C[ns[end]]
+  return temp_O[]
+end
+
+function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteITensorSum)
+  return [expect(ψ, h[(j, j + 1)]) for j in 1:nsites(ψ)]
+end
+
+Sz = [expect(ψ_0, "Sz", n) for n in 1:N]
+energy_infinite = expect(ψ_0, H)
+
+Nfinite = 100
+sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
+Hfinite = MPO(model, sfinite; J=J, h=h)
+ψfinite = randomMPS(sfinite, initstate; linkdims=10)
+@show flux(ψfinite)
+sweeps = Sweeps(15)
+setmaxdim!(sweeps, maxdim)
+setcutoff!(sweeps, cutoff)
+energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps)
+
+nfinite = Nfinite ÷ 2
+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₂),
+)
+
+@show (Sz, Sz_finite)

src/ITensorInfiniteMPS.jl

@@ -35,6 +35,7 @@ include("infinitecanonicalmps.jl")
 include("transfermatrix.jl")
 include("models/models.jl")
 include("models/ising.jl")
+include("models/ising_extended.jl")
 include("models/heisenberg.jl")
 include("models/hubbard.jl")
 include("models/xx.jl")

src/celledvectors.jl

@@ -11,6 +11,7 @@ celltagprefix() = "c="
 celltags(n::Integer) = TagSet(celltagprefix() * string(n))
 celltags(n1::Integer, n2::Integer) = TagSet(celltagprefix() * n1 * "|" * n2)
 
+indextagprefix() = "n="
 #
 # translatecell
 #
@@ -27,6 +28,11 @@ function getcell(ts::TagSet)
   return parse(Int, celltag[(length(celltagprefix()) + 1):end])
 end
 
+function getsite(ts::TagSet)
+  celltag = ITensorInfiniteMPS.tag_starting_with(ts, indextagprefix())
+  return parse(Int, celltag[(length(indextagprefix()) + 1):end])
+end
+
 function translatecell(ts::TagSet, n::Integer)
   ncell = getcell(ts)
   return replacetags(ts, celltags(ncell) => celltags(ncell + n))

src/infinitecanonicalmps.jl

@@ -55,23 +55,23 @@ end
 
 function insert_linkinds!(A; left_dir=ITensors.Out)
   # TODO: use `celllength` here
-  N = length(A)
+  N = nsites(A)
   l = CelledVector{indtype(A)}(undef, N)
   n = N
   s = siteind(A, 1)
   dim = if hasqns(s)
     kwargs = (; dir=left_dir)
     qn_ln = zero_qn(s)
-    [qn_ln => 1]
+    [qn_ln => 1] #Default to 0 on the right
   else
-    kwargs = (;)
+    kwargs = ()
     1
   end
   l[N] = Index(dim, default_link_tags("l", n, 1); kwargs...)
   for n in 1:(N - 1)
     # TODO: is this correct?
     dim = if hasqns(s)
-      qn_ln = (flux(A[n]) + qn_ln) * left_dir
+      qn_ln = flux(A[n]) * left_dir + qn_ln#Fixed a bug on flux conservation
       [qn_ln => 1]
     else
       1
@@ -81,6 +81,9 @@ function insert_linkinds!(A; left_dir=ITensors.Out)
   for n in 1:N
     A[n] = A[n] * onehot(l[n - 1] => 1) * onehot(dag(l[n]) => 1)
   end
+
+  @assert all(i -> flux(i) == zero_qn(s), A) "Flux not invariant under one unit cell translation, not implemented"
+
   return A
 end
 

src/infinitemps.jl

@@ -57,6 +57,20 @@ end
 
 ITensors.siteinds(f::typeof(only), ψ::InfiniteCanonicalMPS) = siteinds(f, ψ.AL)
 
+ITensorInfiniteMPS.getcell(i::Index) = ITensorInfiniteMPS.getcell(tags(i))
+getsite(i::Index) = getsite(tags(i))
+function ITensors.findfirstsiteind(ψ::InfiniteMPS, i::Index)
+  c = ITensorInfiniteMPS.getcell(i)
+  i1 = translatecell(i, -(c - 1))
+  n1 = findfirst(hascommoninds(i1), ψ[Cell(1)])
+  return (c - 1) * nsites(ψ) + n1
+end
+function ITensors.findsites(ψ::InfiniteMPS, T::ITensor)
+  s = filterinds(T; plev=0)
+  return sort([ITensors.findfirstsiteind(ψ, i) for i in s])
+end
+ITensors.findsites(ψ::InfiniteCanonicalMPS, T::ITensor) = findsites(ψ.AL, T)
+
 # For now, only represents nearest neighbor interactions
 # on a linear chain
 struct InfiniteITensorSum
@@ -69,7 +83,29 @@ function Base.getindex(l::InfiniteITensorSum, n1n2::Tuple{Int,Int})
   @assert n2 == n1 + 1
   return l.data[n1]
 end
-nsites(h::InfiniteITensorSum) = length(l.data)
+nsites(h::InfiniteITensorSum) = length(h.data)
+#Gives the range of the Hamiltonian. Useful for better optimized contraction in VUMPS
+nsites_support(h::InfiniteITensorSum) = order.(h.data) ÷ 2
+nsites_support(h::InfiniteITensorSum, n::Int64) = order(h.data[n]) ÷ 2
+
+nrange(h::InfiniteITensorSum) = nrange.(h.data, ncell=nsites(h))
+nrange(h::InfiniteITensorSum, n::Int64) = nrange(h.data[n]; ncell=nsites(h))
+function nrange(h::ITensor; ncell=0)
+  ns = findsites(h; ncell=ncell)
+  return ns[end] - ns[1] + 1
+end
+
+function ITensors.findfirstsiteind(h::ITensor, i::Index, ncell::Int64)
+  c = ITensorInfiniteMPS.getcell(i)
+  n1 = getsite(i)
+  return (c - 1) * ncell + n1
+end
+function ITensors.findsites(h::ITensor; ncell::Int64=1)
+  s = filterinds(h; plev=0)
+  return sort([ITensors.findfirstsiteind(h, i, ncell) for i in s])
+end
+ITensors.findsites(h::InfiniteITensorSum) = [findsites(h, n) for n in 1:nsites(h)]
+ITensors.findsites(h::InfiniteITensorSum, n::Int64) = findsites(h.data[n]; ncell=nsites(h))
 
 ## HDF5 support for the InfiniteCanonicalMPS type
 

src/models/ising.jl

@@ -11,6 +11,19 @@ function ITensors.MPO(::Model"ising", s; J, h)
   return MPO(a, s)
 end
 
+# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ
+function ITensors.OpSum(::Model"ising", n1, n2; J, h)
+  opsum = OpSum()
+  if J != 0
+    opsum += -J, "X", n1, "X", n2
+  end
+  if h != 0
+    opsum += -h / 2, "Z", n1
+    opsum += -h / 2, "Z", n2
+  end
+  return opsum
+end
+
 # H = -J X₁X₂ - h Z₁
 # XXX: use `op` instead of `ITensor`
 function ITensors.ITensor(::Model"ising", s1::Index, s2::Index; J, h)

src/models/ising_extended.jl

@@ -0,0 +1,48 @@
+nrange(::Model"ising_extended") = 3
+# 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 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)
+  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]))
+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