Fix parallel VUMPS and complex number support (#37)

Author: ebertok

src/nullspace.jl

@@ -59,6 +59,7 @@ function blocksparsetensor(blocks::Dict{B,TB}) where {B,TB}
   b1, Tb1 = first(pairs(blocks))
   N = length(b1)
   indstypes = typeof.(inds(Tb1))
+  blocktype = eltype(Tb1)
   indsT = getindex.(indstypes)
   # Determine the indices from the blocks
   for (b, Tb) in pairs(blocks)
@@ -72,7 +73,7 @@ function blocksparsetensor(blocks::Dict{B,TB}) where {B,TB}
       indsTn[bn] = indsTb[n]
     end
   end
-  T = BlockSparseTensor(indsT)
+  T = BlockSparseTensor(blocktype, indsT)
   for (b, Tb) in pairs(blocks)
     if !isempty(Tb)
       T[b] = Tb

src/vumps_localham.jl

@@ -380,10 +380,23 @@ function tdvp_iteration_parallel(
   end
 
   # Sum the Hamiltonian terms in the unit cell
-  for n in 2:Nsites
-    hᴸ[n] = hᴸ[n - 1] * ψ.AL[n] * ψ̃.AL[n] + hᴸ[n]
+  function left_environment_cell(ψ, ψ̃, hᴸ, n)
+    Nsites = nsites(ψ)
+    𝕙ᴸ = copy(hᴸ)
+    for k in reverse((n - Nsites + 2):n)
+      𝕙ᴸ[k] = 𝕙ᴸ[k - 1] * ψ.AL[k] * ψ̃.AL[k] + 𝕙ᴸ[k]
+    end
+    return 𝕙ᴸ[n]
   end
-  Hᴸ = left_environment(hᴸ, ψ; tol=_solver_tol)
+
+  #for k in 2:Nsites
+  #  hᴸ[k] = hᴸ[k - 1] * ψ.AL[k] * ψ̃.AL[k] + hᴸ[k]
+  #end
+  𝕙ᴸ = copy(hᴸ)
+  for k in 1:Nsites
+    𝕙ᴸ[k] = left_environment_cell(ψ, ψ̃, hᴸ, k)
+  end
+  Hᴸ = left_environment(hᴸ, 𝕙ᴸ, ψ; tol=_solver_tol)
 
   for n in 2:Nsites
     hᴿ[n] = hᴿ[n + 1] * ψ.AR[n + 1] * ψ̃.AR[n + 1] + hᴿ[n]
@@ -398,21 +411,27 @@ function tdvp_iteration_parallel(
       Hᴬᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.AL[n] * ψ.C[n], _solver_tol
     )
 
-    C̃[n] = Cvecsₙ[1]
-    Ãᶜ[n] = Avecsₙ[1]
+    C̃[n] = Cvecsₙ
+    Ãᶜ[n] = Avecsₙ
+  end
+
+  function ortho_overlap(AC, C)
+    AL, _ = polar(AC * dag(C), uniqueinds(AC, C))
+    return noprime(AL)
   end
 
-  # TODO: based on minimum singular values of C̃, use more accurate
-  # method for finding Ãᴸ, Ãᴿ
+  function ortho_polar(AC, C)
+    UAC, _ = polar(AC, uniqueinds(AC, C))
+    UC, _ = polar(C, commoninds(C, AC))
+    return noprime(UAC) * noprime(dag(UC))
+  end
+
+
   Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
   Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
   for n in 1:Nsites
-    Ãᴸⁿ, X = polar(Ãᶜ[n] * dag(C̃[n]), uniqueinds(Ãᶜ[n], C̃[n]))
-    Ãᴿⁿ, _ = polar(Ãᶜ[n] * dag(C̃[n - 1]), uniqueinds(Ãᶜ[n], C̃[n - 1]))
-    Ãᴸⁿ = noprime(Ãᴸⁿ)
-    Ãᴿⁿ = noprime(Ãᴿⁿ)
-    Ãᴸ[n] = Ãᴸⁿ
-    Ãᴿ[n] = Ãᴿⁿ
+    Ãᴸ[n] = ortho_polar(Ãᶜ[n], C̃[n])
+    Ãᴿ[n] = ortho_polar(Ãᶜ[n], C̃[n-1])
   end
 
   for n in 1:Nsites