Optimize subspace expansion for long range interactions (#67)

Author: mtfishman

examples/vumps/src/vumps_subspace_expansion.jl

@@ -6,14 +6,16 @@ using ITensorInfiniteMPS
 function tdvp_subspace_expansion(
   H, ψ; time_step, outer_iters, subspace_expansion_kwargs, vumps_kwargs
 )
-  @time for _ in 1:outer_iters
-    println("\nIncrease bond dimension from $(maxlinkdim(ψ))")
+  @time for outer_iter in 1:outer_iters
+    println(
+      "\nIncrease bond dimension $(outer_iter) out of $(outer_iters), starting from dimension $(maxlinkdim(ψ))",
+    )
     println(
       "cutoff = $(subspace_expansion_kwargs[:cutoff]), maxdim = $(subspace_expansion_kwargs[:maxdim])",
     )
-    ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...)
+    ψ = @time subspace_expansion(ψ, H; subspace_expansion_kwargs...)
     println("\nRun VUMPS with new bond dimension $(maxlinkdim(ψ))")
-    ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...)
+    ψ = @time tdvp(H, ψ; time_step=time_step, vumps_kwargs...)
   end
   return ψ
 end

examples/vumps/vumps_2d_heisenberg.jl

@@ -0,0 +1,74 @@
+using ITensorInfiniteMPS
+using ITensors
+
+include(
+  joinpath(
+    pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl"
+  ),
+)
+
+##############################################################################
+# VUMPS parameters
+#
+
+maxdim = 256 # 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-4
+conserve_qns = true
+solver_tol = (x -> x / 10)
+outer_iters = 10 # Number of times to increase the bond dimension
+width = 4
+
+##############################################################################
+# CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
+#
+
+N = width # Number of sites in the unit cell
+
+initstate(n) = isodd(n) ? "↑" : "↓"
+s = infsiteinds("S=1/2", N; conserve_qns, initstate)
+ψ = InfMPS(s, initstate)
+
+function ITensorInfiniteMPS.unit_cell_terms(::Model"heisenberg2D"; width)
+  opsum = OpSum()
+  for i in 1:width
+    # Vertical
+    opsum += -0.5, "S+", i, "S-", mod(i + 1, width)
+    opsum += -0.5, "S-", i, "S+", mod(i + 1, width)
+    opsum += "Sz", i, "Sz", mod(i + 1, width)
+    # Horizontal
+    opsum += -0.5, "S+", i, "S-", i + width
+    opsum += -0.5, "S-", i, "S+", i + width
+    opsum += "Sz", i, "Sz", i + width
+  end
+  return opsum
+end
+model = Model("heisenberg2D")
+
+# Form the Hamiltonian
+H = InfiniteSum{MPO}(model, s; width)
+
+# Check translational invariance
+# println("\nCheck translation invariance of the initial VUMPS state")
+# @show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))
+
+vumps_kwargs = (tol=vumps_tol, maxiter=max_vumps_iters, solver_tol)
+subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
+
+ψ = vumps_subspace_expansion(H, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs)
+
+# Check translational invariance
+# println()
+# println("==============================================================")
+# println()
+# println("\nCheck translation invariance of the final VUMPS state")
+# @show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))
+
+# Sz = [expect(ψ, "Sz", n) for n in 1:N]
+
+energy_infinite = expect(ψ, H)
+@show energy_infinite
+
+## using JLD2
+## jldsave("infmps.jld2"; ψ)

src/subspace_expansion.jl

@@ -26,11 +26,11 @@ function generate_twobody_nullspace(
   end
 
   if range_H == 2
-    ψH2 = noprime(ψ.AL[n1] * H[n1][1] * H[n1][2] * ψ.C[n1] * ψ.AR[n2])
+    ψH2 = (ψ.AL[n1] * ψ.C[n1] * H[n1][1]) * (ψ.AR[n2] * H[n1][2])
+    ψH2 = noprime(ψH2)
   else   # Should be a better version now
     ψH2 =
-      H[n1][end] *
-      (ψ.AR[n2 + range_H - 2] * (ψ′.AR[n2 + range_H - 2] * δʳ(n2 + range_H - 2)))
+      ψ.AR[n2 + range_H - 2] * H[n1][end] * (ψ′.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))
@@ -56,38 +56,40 @@ function generate_twobody_nullspace(
 
     ψH2 = noprime(ψH2)
     for n in 1:(range_H - 2)
-      temp_H2 = δʳ(n2 + range_H - 2 - n)
+      temp_H2_right = δʳ(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] * H[n1 - n][idx] * ψ′.AR[j]
+          temp_H2_right = temp_H2_right * ψ.AR[j] * H[n1 - n][idx] * ψ′.AR[j]
           idx -= 1
         else
-          temp_H2 = temp_H2 * ψ.AR[j] * (δˢ(j) * ψ′.AR[j])
+          temp_H2_right = temp_H2_right * ψ.AR[j] * (δˢ(j) * ψ′.AR[j])
         end
       end
       if common_sites[idx] == n2
-        temp_H2 = temp_H2 * ψ.AR[n2] * H[n1 - n][idx]
+        temp_H2_right = temp_H2_right * ψ.AR[n2] * H[n1 - n][idx]
         idx -= 1
       else
-        temp_H2 = temp_H2 * ψ.AR[n2] * δˢ(n2)
+        temp_H2_right = temp_H2_right * ψ.AR[n2] * δˢ(n2)
       end
       if common_sites[idx] == n1
-        temp_H2 = temp_H2 * (ψ.AL[n1] * ψ.C[n1]) * H[n1 - n][idx]
+        temp_H2_right = temp_H2_right * (ψ.AL[n1] * ψ.C[n1]) * H[n1 - n][idx]
         idx -= 1
       else
-        temp_H2 = temp_H2 * ((ψ.AL[n1] * δˢ(n1)) * ψ.C[n1])
+        temp_H2_right = temp_H2_right * ((ψ.AL[n1] * δˢ(n1)) * ψ.C[n1])
       end
-      for j in 1:n
+      idx = n - idx + 1
+      temp_H2_left = δˡ(n1 - n - 1)
+      for j in reverse(1:n)
         if n1 - j == common_sites[idx]
-          temp_H2 = temp_H2 * ψ.AL[n1 - j] * H[n1 - n][idx] * ψ′.AL[n1 - j]
-          idx -= 1
+          temp_H2_left = temp_H2_left * ψ.AL[n1 - j] * H[n1 - n][idx] * ψ′.AL[n1 - j]
+          idx += 1
         else
-          temp_H2 = temp_H2 * (ψ.AL[n1 - j] * δˢ(n1 - j)) * ψ′.AL[n1 - j]
+          temp_H2_left = temp_H2_left * (ψ.AL[n1 - j] * δˢ(n1 - j)) * ψ′.AL[n1 - j]
         end
       end
-      ψH2 = ψH2 + noprime(temp_H2 * δˡ(n1 - n - 1))
+      ψH2 = ψH2 + noprime(temp_H2_left * temp_H2_right)
     end
   end
   return ψH2

src/vumps_generic.jl

@@ -172,7 +172,7 @@ function tdvp(
     flush(stderr)
   end
   for iter in 1:maxiter
-    ψ, (eᴸ, eᴿ) = tdvp_iteration(
+    iteration_time = @elapsed ψ, (eᴸ, eᴿ) = tdvp_iteration(
       solver,
       ∑h,
       ψ;
@@ -187,7 +187,7 @@ function tdvp(
     maxdimψ = maxlinkdim(ψ[0:(N + 1)])
     if outputlevel > 0
       println(
-        "VUMPS iteration $iter (out of maximum $maxiter). Bond dimension = $maxdimψ, energy = $((eᴸ, eᴿ)), ϵᵖʳᵉˢ = $ϵᵖʳᵉˢ, tol = $tol",
+        "VUMPS iteration $iter (out of maximum $maxiter). Bond dimension = $maxdimψ, energy = $((eᴸ, eᴿ)), ϵᵖʳᵉˢ = $ϵᵖʳᵉˢ, tol = $tol, iteration time = $iteration_time seconds",
       )
       flush(stdout)
       flush(stderr)