More VUMPS optimizations (#65)

* Reorder contraction sequence in the rest of the MPOs in HAC and HC, contracting towards AC

* Replace delta contractions with replaceinds in VUMPS

Author: mtfishman

src/vumps_nonlocalham.jl

@@ -14,62 +14,88 @@ function Base.:*(H::Hᶜ{MPO}, v::ITensor)
   r′ = linkinds(only, ψ′.AR)
   s = siteinds(only, ψ)
   s′ = siteinds(only, ψ′)
-  δˢ(n) = δ(dag(s[n]), prime(s[n]))
-  δʳ(n) = δ(dag(r[n]), prime(r[n]))
-  δˡ(n) = δ(l[n], l′[n])
+
+  # Remove temporarily in favor of `replaceinds`,
+  # which is more efficient.
+  # δˢ(n) = δ(dag(s[n]), prime(s[n]))
+  # δʳ(n) = δ(dag(r[n]), prime(r[n]))
+  # δˡ(n) = δ(l[n], l′[n])
 
   #Build the contribution of the left environment
-  Hᶜᴸv = v * Hᴸ[n] * δʳ(n)
+  # Hᶜᴸv = v * Hᴸ[n] * δʳ(n)
+  Hᶜᴸv = replaceinds(v * Hᴸ[n], r[n] => r[n]')
   #Build the contribution of the right environment
-  Hᶜᴿv = v * δˡ(n) * Hᴿ[n]
+  # Hᶜᴿv = v * δˡ(n) * Hᴿ[n]
+  Hᶜᴿv = replaceinds(v, l[n] => l[n]') * Hᴿ[n]
   #We now start building terms where C overlap with the local Hamiltonian
   # We start with the tensor AL[n] - v - AR[n+1] ... AR[n + range_∑h - 1]
+  # Hᶜʰv =
+  #   δʳ(n + range_∑h - 1) * ψ.AR[n + range_∑h - 1] * ∑h[n][end] * ψ′.AR[n + range_∑h - 1]
   Hᶜʰv =
-    δʳ(n + range_∑h - 1) * ψ.AR[n + range_∑h - 1] * ∑h[n][end] * ψ′.AR[n + range_∑h - 1]
+    replaceinds(ψ.AR[n + range_∑h - 1], r[n + range_∑h - 1] => r[n + range_∑h - 1]') *
+    ∑h[n][end] *
+    ψ′.AR[n + range_∑h - 1]
   common_sites = findsites(ψ, ∑h[n])
   idx = length(∑h[n]) - 1 #list the sites Σh, we start at 2 because n is already taken into account
   for k in reverse(1:(range_∑h - 2))
     if n + k == common_sites[idx]
       Hᶜʰv = Hᶜʰv * ψ.AR[n + k] * ∑h[n][idx] * ψ′.AR[n + k]
       idx -= 1
     else
-      Hᶜʰv = Hᶜʰv * ψ.AR[n + k] * δˢ(n + k) * ψ′.AR[n + k]
+      # Hᶜʰv = Hᶜʰv * ψ.AR[n + k] * δˢ(n + k) * ψ′.AR[n + k]
+      Hᶜʰv = replaceinds(Hᶜʰv * ψ.AR[n + k], s[n + k] => s[n + k]') * ψ′.AR[n + k]
     end
   end
-  Hᶜʰv = v * ψ.AL[n] * δˡ(n - 1) * ∑h[n][1] * ψ′.AL[n] * Hᶜʰv #left extremity
+  # Hᶜʰv = v * ψ.AL[n] * δˡ(n - 1) * ∑h[n][1] * ψ′.AL[n] * Hᶜʰv #left extremity
+  Hᶜʰv = replaceinds(v * ψ.AL[n], l[n - 1] => l[n - 1]') * ∑h[n][1] * ψ′.AL[n] * Hᶜʰv #left extremity
   #Now we are building contributions of the form AL[n - j] ... AL[n] - v - AR[n + 1] ... AR[n + range_∑h - 1 - j]
   for j in 1:(range_∑h - 2)
-    temp_Hᶜʰv = ψ.AL[n - j] * δˡ(n - 1 - j) * ∑h[n - j][1] * ψ′.AL[n - j]
+    # temp_Hᶜʰv = ψ.AL[n - j] * δˡ(n - 1 - j) * ∑h[n - j][1] * ψ′.AL[n - j]
+    temp_Hᶜʰv =
+      replaceinds(ψ.AL[n - j], l[n - 1 - j] => l[n - 1 - j]') * ∑h[n - j][1] * ψ′.AL[n - j]
     common_sites = findsites(ψ, ∑h[n - j])
     idx = 2
     for k in 1:j
       if n - j + k == common_sites[idx]
         temp_Hᶜʰv = temp_Hᶜʰv * ψ.AL[n - j + k] * ∑h[n - j][idx] * ψ′.AL[n - j + k]
         idx += 1
       else
-        temp_Hᶜʰv = temp_Hᶜʰv * ψ.AL[n - j + k] * δˢ(n - j + k) * ψ′.AL[n - j + k]
+        # temp_Hᶜʰv = temp_Hᶜʰv * ψ.AL[n - j + k] * δˢ(n - j + k) * ψ′.AL[n - j + k]
+        temp_Hᶜʰv =
+          replaceinds(temp_Hᶜʰv * ψ.AL[n - j + k], s[n - j + k] => s[n - j + k]') *
+          ψ′.AL[n - j + k]
       end
     end
     # Finished the AL part
-    temp_Hᶜʰv = temp_Hᶜʰv * v
-    for k in (j + 1):(range_∑h - 2)
+    # temp_Hᶜʰv_r =
+    #   ψ.AR[n - j + range_∑h - 1] *
+    #   δʳ(n - j + range_∑h - 1) *
+    #   ∑h[n - j][end] *
+    #   ψ′.AR[n - j + range_∑h - 1]
+    temp_Hᶜʰv_r =
+      replaceinds(
+        ψ.AR[n - j + range_∑h - 1], r[n - j + range_∑h - 1] => r[n - j + range_∑h - 1]'
+      ) *
+      ∑h[n - j][end] *
+      ψ′.AR[n - j + range_∑h - 1]
+    idx = length(∑h[n]) - 1
+    for k in reverse((j + 1):(range_∑h - 2))
       if n - j + k == common_sites[idx]
-        temp_Hᶜʰv = temp_Hᶜʰv * ψ.AR[n - j + k] * ∑h[n - j][idx] * ψ′.AR[n - j + k]
-        idx += 1
+        temp_Hᶜʰv_r = temp_Hᶜʰv_r * ψ.AR[n - j + k] * ∑h[n - j][idx] * ψ′.AR[n - j + k]
+        idx -= 1
       else
-        temp_Hᶜʰv = temp_Hᶜʰv * ψ.AR[n - j + k] * δˢ(n - j + k) * ψ′.AR[n - j + k]
+        # temp_Hᶜʰv_r = temp_Hᶜʰv_r * ψ.AR[n - j + k] * δˢ(n - j + k) * ψ′.AR[n - j + k]
+        temp_Hᶜʰv_r =
+          replaceinds(temp_Hᶜʰv_r * ψ.AR[n - j + k], s[n - j + k] => s[n - j + k]') *
+          ψ′.AR[n - j + k]
       end
     end
-    temp_Hᶜʰv =
-      temp_Hᶜʰv *
-      ψ.AR[n - j + range_∑h - 1] *
-      δʳ(n - j + range_∑h - 1) *
-      ∑h[n - j][end] *
-      ψ′.AR[n - j + range_∑h - 1]
+    temp_Hᶜʰv = temp_Hᶜʰv * v * temp_Hᶜʰv_r
     Hᶜʰv = Hᶜʰv + temp_Hᶜʰv
   end
   Hᶜv = Hᶜᴸv + Hᶜʰv + Hᶜᴿv
-  return Hᶜv * dag(δˡ(n)) * dag(δʳ(n))
+  # return Hᶜv * dag(δˡ(n)) * dag(δʳ(n))
+  return replaceinds(Hᶜv, (l[n], r[n]) => (l[n]', r[n]'))
 end
 
 function Base.:*(H::Hᴬᶜ{MPO}, v::ITensor)
@@ -88,71 +114,101 @@ function Base.:*(H::Hᴬᶜ{MPO}, v::ITensor)
   s = siteinds(only, ψ)
   s′ = siteinds(only, ψ′)
 
-  δˢ(n) = δ(dag(s[n]), prime(s[n]))
-  δʳ(n) = δ(dag(r[n]), prime(r[n]))
-  δˡ(n) = δ(l[n], l′[n])
+  # Remove temporarily in favor of `replaceinds`,
+  # which is more efficient.
+  # δˢ(n) = δ(dag(s[n]), prime(s[n]))
+  # δʳ(n) = δ(dag(r[n]), prime(r[n]))
+  # δˡ(n) = δ(l[n], l′[n])
 
   #Contribution of the left environment
-  Hᴬᶜᴸv = v * Hᴸ[n - 1] * δˢ(n) * δʳ(n)
+  # Hᴬᶜᴸv = v * Hᴸ[n - 1] * δˢ(n) * δʳ(n)
+  Hᴬᶜᴸv = replaceinds(v * Hᴸ[n - 1], (s[n], r[n]) => (s[n]', r[n]'))
   #Contribution of the right environment
-  Hᴬᶜᴿv = v * δˡ(n - 1) * δˢ(n) * Hᴿ[n]
+  # Hᴬᶜᴿv = v * δˡ(n - 1) * δˢ(n) * Hᴿ[n]
+  Hᴬᶜᴿv = replaceinds(v, (l[n - 1], s[n]) => (l[n - 1]', s[n]')) * Hᴿ[n]
   #We now start building terms where AC overlap with the local Hamiltonian
   # We start with the tensor v - AR[n+1] ... AR[n + range_∑h - 1]
+  # Hᴬᶜʰv =
+  #   δʳ(n + range_∑h - 1) * ψ.AR[n + range_∑h - 1] * ∑h[n][end] * ψ′.AR[n + range_∑h - 1] #rightmost extremity
   Hᴬᶜʰv =
-    δʳ(n + range_∑h - 1) * ψ.AR[n + range_∑h - 1] * ∑h[n][end] * ψ′.AR[n + range_∑h - 1] #rightmost extremity
+    replaceinds(ψ.AR[n + range_∑h - 1], r[n + range_∑h - 1] => r[n + range_∑h - 1]') *
+    ∑h[n][end] *
+    ψ′.AR[n + range_∑h - 1] #rightmost extremity
   common_sites = findsites(ψ, ∑h[n])
   idx = length(∑h[n]) - 1  #list the sites Σh, we start at 2 because n is already taken into account
   for k in reverse(1:(range_∑h - 2))
     if n + k == common_sites[idx]
       Hᴬᶜʰv = Hᴬᶜʰv * ψ.AR[n + k] * ∑h[n][idx] * ψ′.AR[n + k]
       idx -= 1
     else
-      Hᴬᶜʰv = Hᴬᶜʰv * ψ.AR[n + k] * δˢ(n + k) * ψ′.AR[n + k]
+      # Hᴬᶜʰv = Hᴬᶜʰv * ψ.AR[n + k] * δˢ(n + k) * ψ′.AR[n + k]
+      Hᴬᶜʰv = replaceinds(Hᴬᶜʰv * ψ.AR[n + k], s[n + k] => s[n + k]') * ψ′.AR[n + k]
     end
   end
-  Hᴬᶜʰv = v * Hᴬᶜʰv * ∑h[n][1] * δˡ(n - 1)
+  # Hᴬᶜʰv = v * Hᴬᶜʰv * ∑h[n][1] * δˡ(n - 1)
+  Hᴬᶜʰv = replaceinds(v * Hᴬᶜʰv * ∑h[n][1], l[n - 1] => l[n - 1]')
   #Now we are building contributions of the form AL[n - j] ... AL[n-1] - v - AR[n + 1] ... AR[n + range_∑h - 1 - j]
   for j in 1:(range_∑h - 1)
-    temp_Hᴬᶜʰv = ψ.AL[n - j] * δˡ(n - j - 1) * ∑h[n - j][1] * ψ′.AL[n - j]
+    # temp_Hᴬᶜʰv = ψ.AL[n - j] * δˡ(n - j - 1) * ∑h[n - j][1] * ψ′.AL[n - j]
+    temp_Hᴬᶜʰv =
+      replaceinds(ψ.AL[n - j], l[n - j - 1] => l[n - j - 1]') * ∑h[n - j][1] * ψ′.AL[n - j]
     common_sites = findsites(ψ, ∑h[n - j])
     idx = 2
     for k in 1:(j - 1)
       if n - j + k == common_sites[idx]
         temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * ψ.AL[n - j + k] * ∑h[n - j][idx] * ψ′.AL[n - j + k]
         idx += 1
       else
-        temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * ψ.AL[n - j + k] * δˢ(n - j + k) * ψ′.AL[n - j + k]
+        # temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * ψ.AL[n - j + k] * δˢ(n - j + k) * ψ′.AL[n - j + k]
+        temp_Hᴬᶜʰv =
+          replaceinds(temp_Hᴬᶜʰv * ψ.AL[n - j + k], s[n - j + k] => s[n - j + k]') *
+          ψ′.AL[n - j + k]
       end
     end
     #Finished with AL, treating the center AC = v
     if j == range_∑h - 1
-      temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * v * ∑h[n - j][end] * δʳ(n - j + range_∑h - 1)
+      # temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * v * ∑h[n - j][end] * δʳ(n - j + range_∑h - 1)
+      temp_Hᴬᶜʰv = replaceinds(
+        temp_Hᴬᶜʰv * v * ∑h[n - j][end], r[n - j + range_∑h - 1] => r[n - j + range_∑h - 1]'
+      )
     else
-      if n == common_sites[idx] #need to check whether we need to branch v
-        temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * v * ∑h[n - j][idx]
-        idx += 1
-      else
-        temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * v * δˢ(n)
-      end
-      for k in (j + 1):(range_∑h - 2)
+      # temp_Hᴬᶜʰv_r =
+      #   ψ.AR[n + range_∑h - 1 - j] *
+      #   δʳ(n - j + range_∑h - 1) *
+      #   ∑h[n - j][end] *
+      #   ψ′.AR[n + range_∑h - 1 - j]
+      temp_Hᴬᶜʰv_r =
+        replaceinds(
+          ψ.AR[n + range_∑h - 1 - j], r[n - j + range_∑h - 1] => r[n - j + range_∑h - 1]'
+        ) *
+        ∑h[n - j][end] *
+        ψ′.AR[n + range_∑h - 1 - j]
+
+      idx = length(∑h[n]) - 1
+      for k in reverse((j + 1):(range_∑h - 2))
         if n + k - j == common_sites[idx]
-          temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * ψ.AR[n + k - j] * ∑h[n - j][idx] * ψ′.AR[n + k - j]
-          idx += 1
+          temp_Hᴬᶜʰv_r = temp_Hᴬᶜʰv_r * ψ.AR[n + k - j] * ∑h[n - j][idx] * ψ′.AR[n + k - j]
+          idx -= 1
         else
-          temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * ψ.AR[n + k - j] * δˢ(n + k - j) * ψ′.AR[n + k - j]
+          # temp_Hᴬᶜʰv_r = temp_Hᴬᶜʰv_r * ψ.AR[n + k - j] * δˢ(n + k - j) * ψ′.AR[n + k - j]
+          temp_Hᴬᶜʰv_r =
+            replaceinds(temp_Hᴬᶜʰv_r * ψ.AR[n + k - j], s[n + k - j] => s[n + k - j]') *
+            ψ′.AR[n + k - j]
         end
       end
-      temp_Hᴬᶜʰv =
-        temp_Hᴬᶜʰv *
-        ψ.AR[n + range_∑h - 1 - j] *
-        δʳ(n - j + range_∑h - 1) *
-        ∑h[n - j][end] *
-        ψ′.AR[n + range_∑h - 1 - j]
+
+      if n == common_sites[idx] #need to check whether we need to branch v
+        temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * v * ∑h[n - j][idx] * temp_Hᴬᶜʰv_r
+      else
+        # temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * v * δˢ(n) * temp_Hᴬᶜʰv_r
+        temp_Hᴬᶜʰv = replaceinds(temp_Hᴬᶜʰv * v * temp_Hᴬᶜʰv_r, s[n] => s[n]')
+      end
     end
     Hᴬᶜʰv = Hᴬᶜʰv + temp_Hᴬᶜʰv
   end
   Hᴬᶜv = Hᴬᶜᴸv + Hᴬᶜʰv + Hᴬᶜᴿv
-  return Hᴬᶜv * dag(δˡ(n - 1)) * dag(δˢ(n)) * dag(δʳ(n))
+  # return Hᴬᶜv * dag(δˡ(n - 1)) * dag(δˢ(n)) * dag(δʳ(n))
+  return replaceinds(Hᴬᶜv, (l[n - 1], s[n], r[n]) => (l[n - 1]', s[n]', r[n]'))
 end
 
 function left_environment(∑h::InfiniteSum{MPO}, ψ::InfiniteCanonicalMPS; tol=1e-15)