General InfiniteBlockMPO working vumps

.gitignore

@@ -2,3 +2,4 @@
 /examples/vumps/Manifest.toml
 /test/Manifest.toml
 *.swp
+.vscode/

Project.toml

@@ -20,10 +20,10 @@ SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 Compat = "3, 4"
 HDF5 = "0.15, 0.16, 0.17"
 ITensorMPS = "0.3"
-ITensors = "0.7, 0.8"
+ITensors = "0.7, 0.8, 0.9"
 Infinities = "0.1"
 IterTools = "1"
-KrylovKit = "0.5, 0.6, 0.7, 0.8, 0.9"
+KrylovKit = "0.5, 0.6, 0.7, 0.8, 0.9, 0.10"
 OffsetArrays = "1"
 QuadGK = "2"
 SplitApplyCombine = "1.2.2"

src/infinitecanonicalmps.jl

@@ -171,8 +171,14 @@ function ITensorMPS.linkinds(ψ::InfiniteMPS)
   return CelledVector([linkinds(ψ, (n, n + 1)) for n in 1:N], translator(ψ))
 end
 
-function InfMPS(s::Vector, f::Function, translator::Function=translatecelltags)
-  return InfMPS(infsiteinds(s, translator), f)
+function InfMPS(
+  eltype::Type{<:Number}, s::Vector, f::Function, translator::Function=translatecelltags; kwargs...
+)
+  return InfMPS(eltype, infsiteinds(s, translator), f; kwargs...)
+end
+
+function InfMPS(s::Vector, f::Function, translator::Function=translatecelltags; kwargs...)
+  return InfMPS(Float64, infsiteinds(s, translator), f; kwargs...)
 end
 
 function indval(iv::Pair)
@@ -185,7 +191,7 @@ function zero_qn(i::Index)
   return zero(qn(first(space(i))))
 end
 
-function insert_linkinds!(A; left_dir=ITensors.Out)
+function insert_linkinds!(A; left_dir=ITensors.Out, extended_linkdim::Int=1)
   # TODO: use `celllength` here
   N = nsites(A)
   l = CelledVector{indtype(A)}(undef, N, translator(A))
@@ -194,19 +200,19 @@ function insert_linkinds!(A; left_dir=ITensors.Out)
   dim = if hasqns(s)
     kwargs = (; dir=left_dir)
     qn_ln = zero_qn(s)
-    [qn_ln => 1] #Default to 0 on the right
+    [qn_ln => extended_linkdim] #Default to 0 on the right
   else
     kwargs = ()
-    1
+    extended_linkdim
   end
   l[N] = Index(dim, default_link_tags("l", n, 1); kwargs...)
   for n in 1:(N - 1)
     # TODO: is this correct?
     dim = if hasqns(s)
       qn_ln = flux(A[n]) * left_dir + qn_ln#Fixed a bug on flux conservation
-      [qn_ln => 1]
+      [qn_ln => extended_linkdim]
     else
-      1
+      extended_linkdim
     end
     l[n] = Index(dim, default_link_tags("l", n, 1); kwargs...)
   end
@@ -219,30 +225,34 @@ function insert_linkinds!(A; left_dir=ITensors.Out)
   return A
 end
 
-function UniformMPS(s::CelledVector, f::Function; left_dir=ITensors.Out)
+function UniformMPS(
+  eltype::Type{<:Number}, s::CelledVector, f::Function; left_dir=ITensors.Out, extended_linkdim::Int=1
+)
   sᶜ¹ = s[Cell(1)]
-  A = InfiniteMPS([ITensor(sⁿ) for sⁿ in sᶜ¹], translator(s))
+  A = InfiniteMPS([ITensor(eltype, sⁿ) for sⁿ in sᶜ¹], translator(s))
   #A.data.translator = translator(s)
   N = length(sᶜ¹)
   for n in 1:N
     Aⁿ = A[n]
     Aⁿ[indval(s[n] => f(n))] = 1.0
     A[n] = Aⁿ
   end
-  insert_linkinds!(A; left_dir=left_dir)
+  insert_linkinds!(A; left_dir=left_dir, extended_linkdim=extended_linkdim)
   return A
 end
 
-function InfMPS(s::CelledVector, f::Function)
+InfMPS(s::CelledVector, f::Function; kwargs...) = InfMPS(Float64, s::CelledVector, f::Function; kwargs...)
+
+function InfMPS(eltype::Type{<:Number}, s::CelledVector, f::Function; extended_linkdim::Int=1)
   # TODO: rename cell_length
   N = length(s)
-  ψL = UniformMPS(s, f; left_dir=ITensors.Out)
-  ψR = UniformMPS(s, f; left_dir=ITensors.In)
+  ψL = UniformMPS(eltype, s, f; left_dir=ITensors.Out, extended_linkdim)
+  ψR = UniformMPS(eltype, s, f; left_dir=ITensors.In, extended_linkdim)
   ψC = InfiniteMPS(N, translator(s))
   l = linkinds(ψL)
   r = linkinds(ψR)
   for n in 1:N
-    ψCₙ = ITensor(dag(l[n])..., r[n]...)
+    ψCₙ = ITensor(eltype, dag(l[n])..., r[n]...)
     ψCₙ[l[n]... => 1, r[n]... => 1] = 1.0
     ψC[n] = ψCₙ
   end

src/vumps_mpo.jl

@@ -103,20 +103,26 @@ function left_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e-
   s = siteinds(only, ψ)
   δʳ(n) = δ(dag(r[n]), prime(r[n]))
   δˡ(n) = δ(l[n], l′[n])
+  # op("Id", s[1]) is permuted w.r.t. below line
   δˢ(n) = δ(dag(s[n]), prime(s[n]))
 
-  eₗ = [0.0]
+  # this code is limited to homogeneous physical spaces on each site within the unit cell!
+  phys_dim = dim(s[1])
+
+  EType = ITensorMPS.promote_itensor_eltype(ψ)
+
+  eₗ = zeros(EType, 1)
   dₕ = size(H[1])[1]
   #Ls = [Vector{ITensor}(undef, dₕ) for j in 1:N]
   Ls = [initialize_left_environment(H, ψ, j; init_last=true) for j in 1:N]
   #Building the L vector for n_1 = 1
   # TM is 2 3 ... N 1
   localR = ψ.C[1] * δʳ(1) * ψ′.C[1] #to revise
-  for n in reverse(1:(dₕ - 1))
+  for a in reverse(1:(dₕ - 1))
     temp_Ls = apply_left_transfer_matrix(
-      translatecell(translator(ψ), Ls[1][n + 1], -1), n + 1, H, ψ, 2 - N
+      translatecell(translator(ψ), Ls[1][a + 1], -1), a + 1, H, ψ, 2 - N
     )
-    for j in 1:n
+    for j in 1:a
       if isassigned(temp_Ls, j)
         if isassigned(Ls[1], j)
           Ls[1][j] += temp_Ls[j]
@@ -125,25 +131,47 @@ function left_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e-
         end
       end
     end
-    if !isempty(H[1][n, n])
-      λ = H[1][n, n][1, 1]
-      δiag = δˢ(1)
+    # starting and ending identity OPs are order 2,
+    # as they have NO link dimension,
+    # because they terminate and do not transport QN numbers.
+    # Intermediate OPs proportional to the identity DO transport QNs and must be order 4!
+    if (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 2)
+      λ = H[1][a, a][1, 1]
+      # δiag = δˢ(1) # not used
       #@assert norm(H[1][n, n] - λ * δiag) == 0 "Non identity diagonal not implemented in MPO"
       @assert abs(λ) <= 1 "Diverging term"
-      eₗ[1] = (Ls[1][n] * localR)[]
-      Ls[1][n] += -(eₗ[1] * denseblocks(δˡ(1)))
+      eₗ[1] = (Ls[1][a] * localR)[]
+      Ls[1][a] += -(eₗ[1] * denseblocks(δˡ(1)))
       if λ == 1
-        A = AOᴸ(ψ, H, n)
-        Ls[1][n], info = linsolve(A, Ls[1][n], 1, -1; tol=tol)
+        A = AOᴸ(ψ, H, a)
+        Ls[1][a], info = linsolve(A, Ls[1][a], 1, -1; tol=tol)
       else
-        println("Not implemented")
+        # println("Not implemented")
+        error(
+          "This case cannot exist! The last and first term of the diagonal of the MPO must be the identity.",
+        )
         flush(stdout)
         flush(stderr)
       end
+    elseif (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 4)
+      # assert diagonal operator!
+      λ = H[1][a, a][1, 1, 1, 1]
+      test_Tensor = reshape(diagm(λ * ones(phys_dim)), (1, phys_dim, phys_dim, 1))
+      # verify that this OP = λ⋅Id
+      @assert norm(test_Tensor - array(H[1][a, a])) ≈ 0.0 "Only operators proportional to the identity are allowed on the diagonal!"
+
+      @assert (λ <= 1.0) && (λ >= 0.0) "The proportionality factor of the operator to the identity must be 0 ≤ λ ≤ 1"
+
+      # solve equations with linsolve
+      eₗ[1] = (Ls[1][a] * localR)[]
+      solo_link = inds(Ls[1][a])[1]
+      Ls[1][a] += -(eₗ[1] * denseblocks(δˡ(1)) * setelt(solo_link[1]))
+      A = AOᴸ(ψ, H, a)
+      Ls[1][a], info = linsolve(A, Ls[1][a], 1, -1; tol=tol)
     end
   end
-  for n in 2:N
-    Ls[n] = apply_local_left_transfer_matrix(Ls[n - 1], H, ψ, n)
+  for a in 2:N
+    Ls[a] = apply_local_left_transfer_matrix(Ls[a - 1], H, ψ, a)
   end
   return CelledVector(Ls), eₗ[1]
 end
@@ -192,7 +220,7 @@ function initialize_right_environment(
     Rs[1] = ITensor(Float64, link, dag(prime(link)))
   end
   for j in 2:(dₕ - 1)
-    mpo_link = only(uniqueinds(H[n - 1][1, j], sit))
+    mpo_link = only(uniqueinds(H[n - 1][dₕ, j], sit))
     Rs[j] = ITensor(Float64, dag(mpo_link), link, dag(prime(link)))
   end
   return Rs
@@ -256,17 +284,22 @@ function right_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e
   δˡ(n) = δ(l[n], dag(prime(l[n])))
   δˢ(n) = δ(dag(s[n]), prime(s[n]))
 
-  eᵣ = [0.0]
+  # this code is limited to homogeneous physical spaces on each site within the unit cell!
+  phys_dim = dim(s[1])
+
+  EType = ITensorMPS.promote_itensor_eltype(ψ)
+
+  eᵣ = zeros(EType, 1)
   dₕ = size(H[1])[1]
   Rs = [initialize_right_environment(H, ψ, j; init_first=true) for j in 1:N]
   #Building the R vector for n_1 = 1
   # TM is 2-N 3-N ... 0
   localL = ψ.C[0] * δˡ(0) * dag(prime(ψ.C[0]))
-  for n in 2:dₕ
+  for a in 2:dₕ
     temp_Rs = apply_right_transfer_matrix(
-      translatecell(translator(ψ), Rs[1][n - 1], 1), n - 1, H, ψ, N
+      translatecell(translator(ψ), Rs[1][a - 1], 1), a - 1, H, ψ, N
     )
-    for j in n:dₕ
+    for j in a:dₕ
       if isassigned(temp_Rs, j)
         if isassigned(Rs[1], j)
           Rs[1][j] += temp_Rs[j]
@@ -275,21 +308,42 @@ function right_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e
         end
       end
     end
-    if !isempty(H[1][n, n])
-      λ = H[1][n, n][1, 1]
-      δiag = δˢ(1)
+    # starting and ending identity OPs are order 2,
+    # as they have NO link dimension,
+    # because they terminate and do not transport QN numbers.
+    # Intermediate OPs proportional to the identity DO transport QNs and must be order 4!
+    if (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 2)
+      λ = H[1][a, a][1, 1]
+      # δiag = δˢ(1)
       #@assert norm(H[1][n, n] - λ * δiag) == 0 "Non identity diagonal not implemented in MPO"
       @assert abs(λ) <= 1 "Diverging term"
-      eᵣ[1] = (localL * Rs[1][n])[]
-      Rs[1][n] += -(eᵣ[1] * denseblocks(δʳ(0)))
+      eᵣ[1] = (localL * Rs[1][a])[]
+      Rs[1][a] += -(eᵣ[1] * denseblocks(δʳ(0)))
       if λ == 1
-        A = AOᴿ(ψ, H, n)
-        Rs[1][n], info = linsolve(A, Rs[1][n], 1, -1; tol=tol)
+        A = AOᴿ(ψ, H, a)
+        Rs[1][a], info = linsolve(A, Rs[1][a], 1, -1; tol=tol)
       else
-        println("Not yet implemented")
+        error(
+          "This case cannot exist! The last and first term of the diagonal of the MPO must be the identity.",
+        )
         flush(stdout)
         flush(stderr)
       end
+    elseif (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 4)
+      # assert diagonal operator!
+      λ = H[1][a, a][1, 1, 1, 1]
+      test_Tensor = reshape(diagm(λ * ones(phys_dim)), (1, phys_dim, phys_dim, 1))
+      # verify that this OP = λ⋅Id
+      @assert norm(test_Tensor - array(H[1][a, a])) ≈ 0.0 "Only operators proportional to the identity are allowed on the diagonal!"
+
+      @assert (λ <= 1.0) && (λ >= 0.0) "The proportionality factor of the operator to the identity must be 0 ≤ λ ≤ 1"
+
+      # solve euqations with linsolve
+      eᵣ[1] = (localL * Rs[1][a])[]
+      solo_link = inds(Rs[1][a])[1]
+      Rs[1][a] += -(eᵣ[1] * denseblocks(δʳ(0)) * setelt(solo_link[1]))
+      A = AOᴿ(ψ, H, a)
+      Rs[1][a], info = linsolve(A, Rs[1][a], 1, -1; tol=tol)
     end
   end
   if N > 1
@@ -365,8 +419,14 @@ function tdvp_iteration_sequential(
   Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, N))
   Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, N))
 
-  eL = zeros(N)
-  eR = zeros(N)
+  EType_ψ = ITensorMPS.promote_itensor_eltype(ψ)
+
+  EType_t = typeof(time_step)
+
+  EType = typeof(one(EType_ψ) * one(EType_t))
+
+  eL = zeros(EType, N)
+  eR = zeros(EType, N)
   for n in 1:N
     L, eL[n] = left_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them
     R, eR[n] = right_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them
@@ -427,8 +487,13 @@ function tdvp_iteration_parallel(
   Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, N))
   Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, N))
 
-  eL = zeros(1)
-  eR = zeros(1)
+  EType_ψ = ITensorMPS.promote_itensor_eltype(ψ)
+  EType_t = typeof(time_step)
+  EType = typeof(one(EType_ψ) * one(EType_t))
+
+  eL = zeros(EType, 1)
+  eR = zeros(EType, 1)
+
   L, eL[1] = left_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them
   R, eR[1] = right_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them
   for n in 1:N