Fix subspace expansion for 1-site unit cell (#44)

Author: LHerviou

src/models/models.jl

@@ -24,7 +24,6 @@ end
 
 function InfiniteITensorSum(model::Model, s::CelledVector; kwargs...)
   N = length(s)
-  H = InfiniteITensorSum(N)
   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)
 end

src/subspace_expansion.jl

@@ -113,61 +113,98 @@ function subspace_expansion(
   NL *= dag(U)
   NR *= dag(V)
 
-  ALⁿ¹, l = ITensors.directsum(
+  ALⁿ¹, newl = ITensors.directsum(
     ψ.AL[n1], dag(NL), uniqueinds(ψ.AL[n1], NL), uniqueinds(NL, ψ.AL[n1]); tags=("Left",)
   )
-  ARⁿ², r = ITensors.directsum(
+  ARⁿ², newr = ITensors.directsum(
     ψ.AR[n2], dag(NR), uniqueinds(ψ.AR[n2], NR), uniqueinds(NR, ψ.AR[n2]); tags=("Right",)
   )
 
-  C = ITensor(dag(l)..., dag(r)...)
+  C = ITensor(dag(newl)..., dag(newr)...)
   ψCⁿ¹ = permute(ψ.C[n1], lⁿ¹..., rⁿ¹...)
   for I in eachindex(ψ.C[n1])
     v = ψCⁿ¹[I]
     if !iszero(v)
       C[I] = ψCⁿ¹[I]
     end
   end
+  if nsites(ψ) == 1
+    ALⁿ² = ITensor(commoninds(C, ARⁿ²), uniqueinds(ALⁿ¹, ψ.AL[n1 - 1])...)
+    il = only(uniqueinds(ALⁿ¹, ALⁿ²))
+    ĩl = only(uniqueinds(ALⁿ², ALⁿ¹))
+    for IV in eachindval(inds(ALⁿ¹))
+      ĨV = replaceind_indval(IV, il => ĩl)
+      v = ALⁿ¹[IV...]
+      if !iszero(v)
+        ALⁿ²[ĨV...] = v
+      end
+    end
+    ARⁿ¹ = ITensor(commoninds(C, ALⁿ¹), uniqueinds(ARⁿ², ψ.AR[n2 + 1])...)
+    ir = only(uniqueinds(ARⁿ², ARⁿ¹))
+    ĩr = only(uniqueinds(ARⁿ¹, ARⁿ²))
+    for IV in eachindval(inds(ARⁿ²))
+      ĨV = replaceind_indval(IV, ir => ĩr)
+      v = ARⁿ²[IV...]
+      if !iszero(v)
+        ARⁿ¹[ĨV...] = v
+      end
+    end
 
-  # Also expand the dimension of the neighboring MPS tensors
-  ALⁿ² = ITensor(dag(l)..., uniqueinds(ψ.AL[n2], ψ.AL[n1])...)
-
-  il = only(uniqueinds(ψ.AL[n2], ALⁿ²))
-  ĩl = only(uniqueinds(ALⁿ², ψ.AL[n2]))
-  for IV in eachindval(inds(ψ.AL[n2]))
-    ĨV = replaceind_indval(IV, il => ĩl)
-    v = ψ.AL[n2][IV...]
-    if !iszero(v)
-      ALⁿ²[ĨV...] = v
+    CL = combiner(newl; tags=tags(only(lⁿ¹)))
+    newind = uniqueinds(CL, newl)
+    newind = replacetags(newind, tags(only(newind)), tags(r[n1 - 1]))
+    temp = δ(dag(newind), newr)
+    ALⁿ² *= CL
+    ALⁿ² *= temp
+    CR = combiner(newr; tags=tags(only(rⁿ¹)))
+    newind = uniqueinds(CR, newr)
+    newind = replacetags(newind, tags(only(newind)), tags(l[n2]))
+    temp = δ(dag(newind), newl)
+    ARⁿ¹ *= CR
+    ARⁿ¹ *= temp
+    C = (C * dag(CL)) * dag(CR)
+    return (ALⁿ¹, ALⁿ²), C, (ARⁿ¹, ARⁿ²)
+  else
+    # Also expand the dimension of the neighboring MPS tensors
+    ALⁿ² = ITensor(dag(newl)..., uniqueinds(ψ.AL[n2], ψ.AL[n1])...)
+
+    il = only(uniqueinds(ψ.AL[n2], ALⁿ²))
+    ĩl = only(uniqueinds(ALⁿ², ψ.AL[n2]))
+    for IV in eachindval(inds(ψ.AL[n2]))
+      ĨV = replaceind_indval(IV, il => ĩl)
+      v = ψ.AL[n2][IV...]
+      if !iszero(v)
+        ALⁿ²[ĨV...] = v
+      end
     end
-  end
 
-  ARⁿ¹ = ITensor(dag(r)..., uniqueinds(ψ.AR[n1], ψ.AR[n2])...)
+    ARⁿ¹ = ITensor(dag(newr)..., uniqueinds(ψ.AR[n1], ψ.AR[n2])...)
 
-  ir = only(uniqueinds(ψ.AR[n1], ARⁿ¹))
-  ĩr = only(uniqueinds(ARⁿ¹, ψ.AR[n1]))
-  for IV in eachindval(inds(ψ.AR[n1]))
-    ĨV = replaceind_indval(IV, ir => ĩr)
-    v = ψ.AR[n1][IV...]
-    if !iszero(v)
-      ARⁿ¹[ĨV...] = v
+    ir = only(uniqueinds(ψ.AR[n1], ARⁿ¹))
+    ĩr = only(uniqueinds(ARⁿ¹, ψ.AR[n1]))
+    for IV in eachindval(inds(ψ.AR[n1]))
+      ĨV = replaceind_indval(IV, ir => ĩr)
+      v = ψ.AR[n1][IV...]
+      if !iszero(v)
+        ARⁿ¹[ĨV...] = v
+      end
     end
-  end
 
-  CL = combiner(l; tags=tags(only(lⁿ¹)))
-  CR = combiner(r; tags=tags(only(rⁿ¹)))
-  ALⁿ¹ *= CL
-  ALⁿ² *= dag(CL)
-  ARⁿ² *= CR
-  ARⁿ¹ *= dag(CR)
-  C = (C * dag(CL)) * dag(CR)
+    CL = combiner(newl; tags=tags(only(lⁿ¹)))
+    CR = combiner(newr; tags=tags(only(rⁿ¹)))
+    ALⁿ¹ *= CL
+    ALⁿ² *= dag(CL)
+    ARⁿ² *= CR
+    ARⁿ¹ *= dag(CR)
+    C = (C * dag(CL)) * dag(CR)
+    return (ALⁿ¹, ALⁿ²), C, (ARⁿ¹, ARⁿ²)
+  end
 
   # TODO: delete or only print when verbose
   ## ψ₂ = ψ.AL[n1] * ψ.C[n1] * ψ.AR[n2]
   ## ψ̃₂ = ALⁿ¹ * C * ARⁿ²
   ## local_energy(ψ, H) = (noprime(ψ * H) * dag(ψ))[]
 
-  return (ALⁿ¹, ALⁿ²), C, (ARⁿ¹, ARⁿ²)
 end
 
 function subspace_expansion(ψ, H; kwargs...)
@@ -181,12 +218,28 @@ function subspace_expansion(ψ, H; kwargs...)
     ALⁿ¹², Cⁿ¹, ARⁿ¹² = subspace_expansion(ψ, H, (n1, n2); kwargs...)
     ALⁿ¹, ALⁿ² = ALⁿ¹²
     ARⁿ¹, ARⁿ² = ARⁿ¹²
-    AL[n1] = ALⁿ¹
-    AL[n2] = ALⁿ²
-    C[n1] = Cⁿ¹
-    AR[n1] = ARⁿ¹
-    AR[n2] = ARⁿ²
+    if N == 1
+      AL[n1] = ALⁿ²
+      C[n1] = Cⁿ¹
+      AR[n2] = ARⁿ¹
+    else
+      AL[n1] = ALⁿ¹
+      AL[n2] = ALⁿ²
+      C[n1] = Cⁿ¹
+      AR[n1] = ARⁿ¹
+      AR[n2] = ARⁿ²
+    end
     ψ = InfiniteCanonicalMPS(AL, C, AR)
   end
   return ψ
 end
+#
+#
+# tt = ψ1.AL[1] * ψ1.C[1]; e = (tt * dag(tt))[]
+# @show norm( ψ1.AL[1] * ψ1.C[1] -  ψ1.C[0] * ψ1.AR[1] )
+# l = linkinds(only, ψ1.AL)
+# r = linkinds(only, ψ1.AR)
+# tt = δ(l[0], prime(dag(l[0]))) * ψ1.AL[1] * dag(prime(ψ1)).AL[1] * dag(δ(s[1], dag(prime(s[1]))));
+# tt == denseblocks(δ(l[1], prime(dag(l[1]))))
+# tt = δ(dag(r[1]), prime(r[1])) * ψ1.AR[1] * dag(prime(ψ1)).AR[1] * dag(δ(s[1], dag(prime(s[1]))));
+# tt == denseblocks(δ(dag(r[0]), prime(r[0])))

test/test_vumps.jl

@@ -38,15 +38,15 @@ using Random
 
   for multisite_update_alg in ["sequential", "parallel"],
     conserve_qns in [true, false],
-    nsites in [2, 3, 4],
+    nsites in [1, 2, 3, 4],
     time_step in [-Inf, -0.5]
 
-    if conserve_qns
-      # Flux density is inconsistent for odd unit cells
-      isodd(nsites) && continue
-    end
+    #if conserve_qns
+    # Flux density is inconsistent for odd unit cells
+    #  isodd(nsites) && continue
+    #end
 
-    space_ = fill(space_shifted(model, 0; conserve_qns=conserve_qns), nsites)
+    space_ = fill(space_shifted(model, 1; conserve_qns=conserve_qns), nsites)
     s = infsiteinds("S=1/2", nsites; space=space_)
     ψ = InfMPS(s, initstate)
 
@@ -123,7 +123,7 @@ end
 ##     space = (("SzParity", 1, 2) => χ ÷ 2) ⊕ (("SzParity", 0, 2) => χ ÷ 2)
 ##     ψ = InfiniteMPS(ComplexF64, s; space=space)
 ##     randn!.(ψ)
-## 
+##
 ##     ψ = orthogonalize(ψ, :)
 ##     @test prod(ψ.AL[1:N]) * ψ.C[N] ≈ ψ.C[0] * prod(ψ.AR[1:N])
 ##   end