Replace the local ITensors by MPOs in InfiniteITensorSum (#45)

Author: LHerviou

examples/vumps/vumps_ising_extended.jl

@@ -44,39 +44,6 @@ for j in 1:outer_iters
   ψ_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)
 

examples/vumps/vumps_ising_noncontiguous.jl

@@ -31,7 +31,11 @@ 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]))
+  temp = MPO(3)
+  temp[1] = H[x][1]
+  temp[2] = H[x][2]
+  temp[3] = delta(prime(s[x + 3]), dag(s[x + 3]))
+  H.data[x] = temp
 end
 
 @show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...));
@@ -47,39 +51,6 @@ for j in 1:outer_iters
   ψ_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)
 

src/ITensors.jl

@@ -72,6 +72,10 @@ function Base.startswith(tag::Tag, subtag::Tag)
 end
 
 function tag_starting_with(ts::TagSet, prefix)
+  x = findfirst(t -> startswith(t, Tag(prefix)), ts)
+  if isnothing(x) #in case the function is called one a link leg. Can be probably improved
+    return x
+  end
   return ts[findfirst(t -> startswith(t, Tag(prefix)), ts)]
 end
 

src/celledvectors.jl

@@ -25,16 +25,22 @@ indextagprefix() = "n="
 # Determine the cell `n` from the tag `"c=n"`
 function getcell(ts::TagSet)
   celltag = tag_starting_with(ts, celltagprefix())
+  if isnothing(celltag) #dealing with link legs
+    return celltag
+  end
   return parse(Int, celltag[(length(celltagprefix()) + 1):end])
 end
 
 function getsite(ts::TagSet)
-  celltag = ITensorInfiniteMPS.tag_starting_with(ts, indextagprefix())
+  celltag = tag_starting_with(ts, indextagprefix())
   return parse(Int, celltag[(length(indextagprefix()) + 1):end])
 end
 
 function translatecell(ts::TagSet, n::Integer)
   ncell = getcell(ts)
+  if isnothing(ncell)
+    return ts
+  end
   return replacetags(ts, celltags(ncell) => celltags(ncell + n))
 end
 
@@ -49,6 +55,7 @@ function translatecell(is::Union{<:Tuple,<:Vector}, n::Integer)
 end
 
 translatecell(T::ITensor, n::Integer) = ITensors.setinds(T, translatecell(inds(T), n))
+translatecell(T::MPO, n::Integer) = translatecell.(T, n)
 
 struct CelledVector{T} <: AbstractVector{T}
   data::Vector{T}

src/infinitecanonicalmps.jl

@@ -115,3 +115,37 @@ function InfMPS(s::CelledVector, f::Function)
   end
   return ψ = InfiniteCanonicalMPS(ψL, ψC, ψR)
 end
+
+function ITensors.expect(ψ::InfiniteCanonicalMPS, o, n)
+  s = siteinds(only, ψ.AL)
+  return (noprime(ψ.AL[n] * ψ.C[n] * op(o, s[n])) * dag(ψ.AL[n] * ψ.C[n]))[]
+end
+
+function ITensors.expect(ψ::InfiniteCanonicalMPS, h::MPO)
+  l = linkinds(ITensorInfiniteMPS.only, ψ.AL)
+  r = linkinds(ITensorInfiniteMPS.only, ψ.AR)
+  s = siteinds(ITensorInfiniteMPS.only, ψ)
+  δˢ(n) = ITensorInfiniteMPS.δ(dag(s[n]), prime(s[n]))
+  δˡ(n) = ITensorInfiniteMPS.δ(l[n], prime(dag(l[n])))
+  δʳ(n) = ITensorInfiniteMPS.δ(dag(r[n]), prime(r[n]))
+  ψ′ = prime(dag(ψ))
+
+  ns = ITensorInfiniteMPS.findsites(ψ, h)
+  nrange = ns[end] - ns[1] + 1
+  idx = 2
+  temp_O = δˡ(ns[1] - 1) * ψ.AL[ns[1]] * ψ′.AL[ns[1]] * h[1]
+  for n in (ns[1] + 1):(ns[1] + nrange - 1)
+    if n == ns[idx]
+      temp_O = temp_O * ψ.AL[n] * ψ′.AL[n] * h[idx]
+      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]) for j in 1:nsites(ψ)]
+end

src/infinitemps.jl

@@ -57,14 +57,24 @@ end
 
 ITensors.siteinds(f::typeof(only), ψ::InfiniteCanonicalMPS) = siteinds(f, ψ.AL)
 
-ITensorInfiniteMPS.getcell(i::Index) = ITensorInfiniteMPS.getcell(tags(i))
+getcell(i::Index) = ITensorInfiniteMPS.getcell(tags(i))
 getsite(i::Index) = getsite(tags(i))
+
+#Before, for a two site Hamiltonian, findsites(ψ, H[3]) would return [1, 2]
+#Appeared unsafe for index purpose
 function ITensors.findfirstsiteind(ψ::InfiniteMPS, i::Index)
   c = ITensorInfiniteMPS.getcell(i)
-  i1 = translatecell(i, -(c - 1))
-  n1 = findfirst(hascommoninds(i1), ψ[Cell(1)])
+  n1 = getsite(i)
   return (c - 1) * nsites(ψ) + n1
 end
+
+function ITensors.findsites(ψ::InfiniteMPS, T::MPO)
+  s = [noprime(filterinds(T[x]; plev=1)[1]) for x in 1:length(T)]
+  return sort([ITensors.findfirstsiteind(ψ, i) for i in s])
+end
+ITensors.findsites(ψ::InfiniteCanonicalMPS, T::MPO) = findsites(ψ.AL, T)
+
+#Kept for historical reason
 function ITensors.findsites(ψ::InfiniteMPS, T::ITensor)
   s = filterinds(T; plev=0)
   return sort([ITensors.findfirstsiteind(ψ, i) for i in s])
@@ -74,38 +84,57 @@ ITensors.findsites(ψ::InfiniteCanonicalMPS, T::ITensor) = findsites(ψ.AL, T)
 # For now, only represents nearest neighbor interactions
 # on a linear chain
 struct InfiniteITensorSum
-  data::CelledVector{ITensor}
+  data::CelledVector{MPO}
 end
 InfiniteITensorSum(N::Int) = InfiniteITensorSum(Vector{ITensor}(undef, N))
-InfiniteITensorSum(data::Vector{ITensor}) = InfiniteITensorSum(CelledVector(data))
+InfiniteITensorSum(data::Vector{MPO}) = InfiniteITensorSum(CelledVector(data))
+InfiniteITensorSum(data::Vector{ITensor}) = InfiniteITensorSum(CelledVector(MPO.(data)))
 function Base.getindex(l::InfiniteITensorSum, n1n2::Tuple{Int,Int})
   n1, n2 = n1n2
   @assert n2 == n1 + 1
   return l.data[n1]
 end
+function Base.getindex(l::InfiniteITensorSum, n1::Int)
+  return l.data[n1]
+end
 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
+nsites_support(h::InfiniteITensorSum) = length.(h.data)
+nsites_support(h::InfiniteITensorSum, n::Int64) = length(h.data[n])
 
 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)
+function nrange(h::MPO; ncell=1)
   ns = findsites(h; ncell=ncell)
   return ns[end] - ns[1] + 1
 end
 
-function ITensors.findfirstsiteind(h::ITensor, i::Index, ncell::Int64)
+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))
+function ITensors.findfirstsiteind(i::Index, ncell::Int64)
   c = ITensorInfiniteMPS.getcell(i)
   n1 = getsite(i)
   return (c - 1) * ncell + n1
 end
+function ITensors.findsites(h::MPO; ncell::Int64=1)
+  s = [filterinds(h[x]; plev=1)[1] for x in 1:length(h)]
+  return sort([ITensors.findfirstsiteind(i, ncell) for i in s])
+end
+
+#Kept for historical reasons
+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 = 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/models.jl

@@ -24,13 +24,18 @@ end
 
 function InfiniteITensorSum(model::Model, s::CelledVector; kwargs...)
   N = length(s)
-  tensors = [ITensor(model, s, n; kwargs...) for n in 1:N] #slightly improved version. Note: the current implementation does not really allow for staggered potentials for example
-  return InfiniteITensorSum(tensors)
+  mpos = [MPO(model, s, n; kwargs...) for n in 1:N] #slightly improved version. Note: the current implementation does not really allow for staggered potentials for example
+  return InfiniteITensorSum(mpos)
 end
 
-# Version accepting IndexSet
-function ITensors.ITensor(model::Model, s::CelledVector, n::Int64; kwargs...)
+# MPO building version
+function ITensors.MPO(model::Model, s::CelledVector, n::Int64; kwargs...)
   n1, n2 = 1, 2
   opsum = OpSum(model, n1, n2; kwargs...)
-  return prod(MPO(opsum, [s[x] for x in n:(n + nrange(model) - 1)])) #modification to allow for more than two sites per term in the Hamiltonians
+  return MPO(opsum, [s[x] for x in n:(n + nrange(model) - 1)]) #modification to allow for more than two sites per term in the Hamiltonians
+end
+
+# Version accepting IndexSet
+function ITensors.ITensor(model::Model, s::CelledVector, n::Int64; kwargs...)
+  return prod(MPO(model, s, n; kwargs...)) #modification to allow for more than two sites per term in the Hamiltonians
 end

src/subspace_expansion.jl

@@ -45,48 +45,61 @@ function subspace_expansion(
   nL = uniqueinds(NL, ψ.AL[n1])
   nR = uniqueinds(NR, ψ.AR[n2])
   if range_H == 2
-    ψH2 = noprime(ψ.AL[n1] * H[(n1, n2)] * ψ.C[n1] * ψ.AR[n2])
+    ψH2 = noprime(ψ.AL[n1] * H[n1][1] * H[n1][2] * ψ.C[n1] * ψ.AR[n2])
   else   # Should be a better version now
     ψH2 =
-      H[(n1, n2)] * ψ.AR[n2 + range_H - 2] * ψ′.AR[n2 + range_H - 2] * δʳ(n2 + range_H - 2)
+      H[n1][end] * ψ.AR[n2 + range_H - 2] * ψ′.AR[n2 + range_H - 2] * δʳ(n2 + range_H - 2)
     common_sites = findsites(ψ, H[(n1, n2)])
     idx = length(common_sites) - 1
     for j in reverse(1:(range_H - 3))
       if n2 + j == common_sites[idx]
-        ψH2 = ψH2 * ψ.AR[n2 + j] * ψ′.AR[n2 + j]
+        ψH2 = ψH2 * ψ.AR[n2 + j] * ψ′.AR[n2 + j] * H[n1][idx]
         idx -= 1
       else
         ψH2 = ψH2 * ψ.AR[n2 + j] * δˢ(n2 + j) * ψ′.AR[n2 + j]
       end
     end
-    ψH2 = noprime(ψH2 * ψ.AL[n1] * ψ.C[n1] * ψ.AR[n2])
+    if common_sites[idx] == n2
+      ψH2 = ψH2 * ψ.AR[n2] * H[n1][idx]
+      idx -= 1
+    else
+      ψH2 = ψH2 * ψ.AR[n2] * δˢ(n2)
+    end
+    if common_sites[idx] == n1
+      ψH2 = ψH2 * ψ.AL[n1] * ψ.C[n1] * H[n1][idx]
+      idx -= 1
+    else
+      ψH2 = ψH2 * ψ.AL[n1] * ψ.C[n1] * δˢ(n1)
+    end
+
+    ψH2 = noprime(ψH2)
     for n in 1:(range_H - 2)
-      temp_H2 = H[(n1 - n, n2 - n)] * δʳ(n2 + range_H - 2 - n)
-      common_sites = findsites(ψ, H[(n1 - n, n2 - n)])
+      temp_H2 = δʳ(n2 + range_H - 2 - n)
+      common_sites = findsites(ψ, H[n1 - n])
       idx = length(common_sites)
       for j in (n2 + range_H - 2 - n):-1:(n2 + 1)
         if j == common_sites[idx]
-          temp_H2 = temp_H2 * ψ.AR[j] * ψ′.AR[j]
+          temp_H2 = temp_H2 * ψ.AR[j] * ψ′.AR[j] * H[n1 - n][idx]
           idx -= 1
         else
           temp_H2 = temp_H2 * ψ.AR[j] * ψ′.AR[j] * δˢ(j)
         end
       end
       if common_sites[idx] == n2
-        temp_H2 = temp_H2 * ψ.AR[n2]
+        temp_H2 = temp_H2 * ψ.AR[n2] * H[n1 - n][idx]
         idx -= 1
       else
         temp_H2 = temp_H2 * ψ.AR[n2] * δˢ(n2)
       end
       if common_sites[idx] == n1
-        temp_H2 = temp_H2 * ψ.AL[n1] * ψ.C[n1]
+        temp_H2 = temp_H2 * ψ.AL[n1] * ψ.C[n1] * H[n1 - n][idx]
         idx -= 1
       else
         temp_H2 = temp_H2 * ψ.AL[n1] * ψ.C[n1] * δˢ(n1)
       end
       for j in 1:n
         if n1 - j == common_sites[idx]
-          temp_H2 = temp_H2 * ψ.AL[n1 - j] * ψ′.AL[n1 - j]
+          temp_H2 = temp_H2 * ψ.AL[n1 - j] * ψ′.AL[n1 - j] * H[n1 - n][idx]
           idx -= 1
         else
           temp_H2 = temp_H2 * ψ.AL[n1 - j] * δˢ(n1 - j) * ψ′.AL[n1 - j]