diff --git a/.gitignore b/.gitignore
index ace59e2..b73f6ce 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,4 @@
 /examples/vumps/Manifest.toml
 /test/Manifest.toml
 *.swp
+.vscode/
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index dce4615..25ef896 100644
--- a/Project.toml
+++ b/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"
diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl
index 6d1dd6a..b3f5196 100644
--- a/src/infinitecanonicalmps.jl
+++ b/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,9 +225,11 @@ 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
@@ -229,20 +237,22 @@ function UniformMPS(s::CelledVector, f::Function; left_dir=ITensors.Out)
     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
diff --git a/src/vumps_mpo.jl b/src/vumps_mpo.jl
index 4a7c561..abc2ff8 100644
--- a/src/vumps_mpo.jl
+++ b/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
diff --git a/test/test_vumps_GeneralMPO.jl b/test/test_vumps_GeneralMPO.jl
new file mode 100644
index 0000000..34d5629
--- /dev/null
+++ b/test/test_vumps_GeneralMPO.jl
@@ -0,0 +1,305 @@
+using ITensors, ITensorMPS
+using ITensorInfiniteMPS
+
+base_path = joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src")
+src_files = ["vumps_subspace_expansion.jl", "entropy.jl"]
+for f in src_files
+  include(joinpath(base_path, f))
+end
+
+"""
+build the transverse field Ising model with exponential interactions in the thermodynamic limit 𝑁 → ∞
+use convention H = -J ∑_{n,m=−∞}^∞ σˣₙ σˣₘ ⋅ λ^(-|n-m-1|)  - hz ∑_{n=−∞}^∞ σᶻₙ
+with 
+"""
+function InfiniteExpHTFI(
+  sites::CelledVector{<:Index},
+  λ, # exponential decay base
+  J, # interaction kinetic (tunneling) coupling
+  hz; # interaction strength of the transverse field
+  kwargs...,
+)
+  if abs(λ) > 1.0
+    throw(
+      ArgumentError("cannot implement exponential decay with base larger than 1, λ = $(λ)!")
+    )
+  end
+
+  link_dimension = 3
+
+  N = length(sites)
+
+  EType = eltype(union(λ, J, hz))
+
+  linkindices = CelledVector(
+    if hasqns(sites.data)
+      [Index([QN("SzParity", 1, 2) => 1], "Link,c=1,n=$n") for n in 1:N]
+    else
+      [Index(1, "Link,c=1,n=$(n)") for n in 1:N]
+    end,
+  )
+
+  mpos = [Matrix{ITensor}(undef, link_dimension, link_dimension) for i in 1:N]
+  for n in 1:N
+    # define local matrix Hmat with empty tensors as local operators
+    Hmat = fill(
+      ITensor(EType, dag(sites[n]), prime(sites[n])), link_dimension, link_dimension
+    )
+    # left link index ll with daggered QN conserving direction (if applicable)
+    ll = dag(linkindices[n - 1])
+    # right link index rl
+    rl = linkindices[n]
+
+    # add both Identities as netral elements in the MPO
+    # replace all known tensors from empty to known interactions
+    Hmat[1, 1] = op("Id", sites[n])
+    Hmat[3, 3] = op("Id", sites[n])
+    # local nearest neighbour and exp. decaying interaction terms
+    Hmat[2, 1] = op("X", sites[n])
+    if !iszero(λ)
+      Hmat[2, 2] = op("Id", sites[n]) * λ  # λ Id,  on the diagonal
+    end
+    Hmat[3, 2] = op("X", sites[n]) * -J # Jxx σˣ
+    if !iszero(hz)
+      Hmat[3, 1] = op("Z", sites[n]) * -hz # hz σᶻ
+    end
+
+    # add all missing links that connect the interaction
+    # operators in the unit cell
+    Hmat[2, 1] = setelt(ll[1]) * Hmat[2, 1]
+    Hmat[1, 2] = Hmat[1, 2] * setelt(rl[1])
+    Hmat[2, 2] = setelt(ll[1]) * Hmat[2, 2] * setelt(rl[1])
+    Hmat[3, 2] = setelt(rl[1]) * Hmat[3, 2]
+    Hmat[2, 3] = setelt(ll[1]) * Hmat[2, 3]
+    mpos[n] = Hmat
+  end
+
+  return InfiniteBlockMPO(mpos, sites.translator)
+end
+
+# sanity check if I can construct the MPO for the normal TFI correctly
+"""
+build the transverse field Ising model with nearest neighbour interactions in the thermodynamic limit 𝑁 → ∞
+use convention H = -J ∑_{n=−∞}^∞ σˣₙ σˣₙ₊₁  - hz ∑_{n=−∞}^∞ σᶻₙ
+with 
+"""
+function InfiniteHTFI(
+  sites::CelledVector{<:Index},
+  kinetic_coupling,
+  hz; # interaction kinetic (tunneling) coupling
+  kwargs...,
+)
+  if abs(λ) > 1.0
+    throw(
+      ArgumentError("cannot implement exponential decay with base larger than 1, λ = $(λ)!")
+    )
+  end
+  link_dimension = 3
+
+  N = length(sites)
+
+  EType = eltype(union(J, hz))
+
+  linkindices = CelledVector(
+    if hasqns(sites.data)
+      [Index([QN("SzParity", 1, 2) => 1], "Link,c=1,n=$n") for n in 1:N]
+    else
+      [Index(1, "Link,c=1,n=$(n)") for n in 1:N]
+    end,
+  )
+
+  mpos = [Matrix{ITensor}(undef, link_dimension, link_dimension) for i in 1:N]
+  for n in 1:N
+    # define local matrix Hmat with empty tensors as local operators
+    Hmat = fill(
+      ITensor(EType, dag(sites[n]), prime(sites[n])), link_dimension, link_dimension
+    )
+    # left link index ll with daggered QN conserving direction (if applicable)
+    ll = dag(linkindices[n - 1])
+    # right link index rl
+    rl = linkindices[n]
+
+    # add both Identities as netral elements in the MPO
+    # replace all known tensors from empty to known interactions
+    Hmat[1, 1] = op("Id", sites[n])
+    Hmat[3, 3] = op("Id", sites[n])
+    # local nearest neighbour and exp. decaying interaction terms
+    Hmat[2, 1] = op("X", sites[n])
+    Hmat[3, 2] = op("X", sites[n]) * -J # Jxx σˣ
+    if !iszero(hz)
+      Hmat[3, 1] = op("Z", sites[n]) * -hz # hz σᶻ
+    end
+
+    # add all missing links that connect the 
+    # interaction operators in the unit cell
+    Hmat[2, 1] = setelt(ll[1]) * Hmat[2, 1]
+    # Hmat[1,2] = Hmat[1,2] * setelt(rl[1])
+    Hmat[3, 2] = setelt(rl[1]) * Hmat[3, 2]
+    # Hmat[2,3] = setelt(ll[1]) * Hmat[2,3]
+    mpos[n] = Hmat
+  end
+  return InfiniteBlockMPO(mpos, sites.translator)
+end
+
+function expect_two_site(ψ::InfiniteCanonicalMPS, h::ITensor, n1n2)
+  n1, n2 = n1n2
+  ϕ = ψ.AL[n1] * ψ.AL[n2] * ψ.C[n2]
+  return inner(ϕ, apply(h, ϕ))
+end
+
+function expect_two_site(ψ::InfiniteCanonicalMPS, h::MPO, n1n2)
+  return expect_two_site(ψ, contract(h), n1n2)
+end
+
+function energy_local(ψ1, ψ2, h::ITensor)
+  ϕ = ψ1 * ψ2
+  return (noprime(ϕ * h) * dag(ϕ))[]
+end
+
+energy_local(ψ1, ψ2, h::MPO) = energy_local(ψ1, ψ2, prod(h))
+# Check translational invariance
+function ITensorMPS.expect(ψ, o)
+  return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[]
+end
+
+maxdim = 16 # Maximum bond dimension
+cutoff = 1e-10 # Singular value cutoff when increasing the bond dimension
+max_vumps_iters = 100 # Maximum number of iterations of the VUMPS/TDVP algorithm at a fixed bond dimension
+tol = 1e-8 # Precision error tolerance for outer loop of VUMPS or TDVP
+outer_iters = 4 # Number of times to increase the bond dimension
+time_step = -Inf # -Inf corresponds to VUMPS, finite time_step corresponds to TDVP
+solver_tol = (x -> x / 100) # Tolerance for the local solver (eigsolve in VUMPS and exponentiate in TDVP)
+multisite_update_alg = "parallel" # Choose between ["sequential", "parallel"]. Only parallel works with TDVP.
+# conserve_qns = true # Whether or not to conserve spin parity
+
+# for conserve_qns in [true, false]
+
+println("Test with Sx conservation = $(conserve_qns)")
+
+nsite = 2 # Number of sites in the unit cell
+initstate(n) = "↑"
+s = infsiteinds("S=1/2", nsite; initstate, conserve_szparity=conserve_qns)
+
+ψ = InfMPS(s, initstate)
+
+@show norm(contract(ψ.AL[1:nsite]..., ψ.C[nsite]) - contract(ψ.C[0], ψ.AR[1:nsite]...))
+
+# J  = -0.25
+# hz = 2.0
+# λ  = 0.4
+
+J = 1
+hz = 1.1
+
+# H_test = InfiniteHTFI(s,J,hz)
+H_test = InfiniteHTFI(s, J, hz)
+H_test0 = InfiniteExpHTFI(s, 0.0, J, hz)
+
+vumps_kwargs = (
+  tol=tol,
+  maxiter=max_vumps_iters,
+  solver_tol=solver_tol,
+  multisite_update_alg=multisite_update_alg,
+)
+
+subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
+
+H_ref = InfiniteSum{MPO}(Model("ising"), s; J=J, h=hz)
+
+ψ_ref = vumps_subspace_expansion(
+  H_ref, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs
+)
+
+ψ_test_NN = vumps_subspace_expansion(
+  H_test, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs
+)
+ψ_test0_NN = vumps_subspace_expansion(
+  H_test0, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs
+)
+
+E_ref = energy_local(ψ_ref.AL[1], ψ_ref.AL[2] * ψ_ref.C[2], H_ref[(1, 2)])
+
+L_test, energy_test = ITensorInfiniteMPS.left_environment(H_test, ψ_test_NN; tol=1e-10);
+L_test0, energy_test0 = ITensorInfiniteMPS.left_environment(H_test0, ψ_test0_NN; tol=1e-10);
+
+energy_test /= length(ψ_test_NN)
+energy_test0 /= length(ψ_test0_NN)
+
+@test energy_test ≈ E_ref
+@test energy_test0 ≈ E_ref
+
+E_exact = reference(Model("ising"), Observable("energy"); h=hz, J)
+
+@test isapprox(E_ref, E_exact; rtol=1e-10)
+@test isapprox(energy_test, E_exact; rtol=1e-10)
+@test isapprox(energy_test0, E_exact; rtol=1e-10)
+
+# energy_local(ψ_NN.AL[1], ψ_NN.AL[2] * ψ.C[2], H_test[(1, 2)])
+# energy_local(ψ_test_NN,ψ_test_NN,H_test)
+
+function ITensorMPS.expect(ψ, o)
+  return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[]
+end
+
+@test isapprox(
+  expect(ψ_ref.AL[1] * ψ_ref.C[1], "Z"),
+  expect(ψ_test_NN.AL[1] * ψ_test_NN.C[1], "Z");
+  atol=1e-7,
+)
+@test isapprox(
+  abs(expect(ψ_ref.AL[1] * ψ_ref.C[1], "X")),
+  abs(expect(ψ_test_NN.AL[1] * ψ_test_NN.C[1], "X"));
+  atol=1e-7,
+)
+
+@test isapprox(
+  expect(ψ_ref.AL[1] * ψ_ref.C[1], "Z"),
+  expect(ψ_test0_NN.AL[1] * ψ_test0_NN.C[1], "Z");
+  atol=1e-7,
+)
+@test isapprox(
+  abs(expect(ψ_ref.AL[1] * ψ_ref.C[1], "X")),
+  abs(expect(ψ_test0_NN.AL[1] * ψ_test0_NN.C[1], "X"));
+  atol=1e-7,
+)
+
+entropy(ψ_ref, 1)
+entropy(ψ_test_NN, 1)
+entropy(ψ_test0_NN, 1)
+
+@test isapprox(entropy(ψ_ref, 1), entropy(ψ_test_NN, 1); rtol=1e-7)
+@test isapprox(entropy(ψ_ref, 1), entropy(ψ_test0_NN, 1); rtol=1e-7)
+
+##################
+#### actually exponential interactions
+#################
+λ = 0.6
+s = infsiteinds("S=1/2", nsite; initstate, conserve_szparity=false)
+
+vumps_kwargs = (
+  tol=1e-6,
+  maxiter=max_vumps_iters,
+  solver_tol=solver_tol,
+  multisite_update_alg=multisite_update_alg,
+)
+ψ = InfMPS(s, initstate)
+
+H_test_exp = InfiniteExpHTFI(s, λ, J, hz)
+
+ψ_exp = vumps_subspace_expansion(
+  H_test_exp, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs
+)
+
+L_test, energy_test_exp_L = ITensorInfiniteMPS.left_environment(H_test_exp, ψ_exp; tol=1e-6);
+R_test, energy_test_exp_R = ITensorInfiniteMPS.right_environment(
+  H_test_exp, ψ_exp; tol=1e-6
+);
+
+expect(ψ_exp)
+
+E_gs_λ04_J1_h1_MPSKit = -1.8177300191088515 # computed at χ=16, λ=0.4, Jₓₓ=1.0, hz=1.0
+E_gs_λ04_J1_h11_MPSKit = -1.8497463013720623 # computed at χ=16, λ=0.4, Jₓₓ=1.0, hz=1.1
+E_gs_04_ITensor = energy_test_exp_L / length(ψ_exp)
+@test E_gs_λ04_J1_h11_MPSKit ≈ E_gs_04_ITensor
+
+# end