diff --git a/src/ITensorInfiniteMPS.jl b/src/ITensorInfiniteMPS.jl
index 1c17184..e147a75 100644
--- a/src/ITensorInfiniteMPS.jl
+++ b/src/ITensorInfiniteMPS.jl
@@ -35,7 +35,7 @@ include("celledvectors.jl")
 include("abstractinfinitemps.jl")
 include("infinitemps.jl")
 include("infinitecanonicalmps.jl")
-include("infinitempomatrix.jl")
+include("infinitempoblock.jl")
 include("infinitempo.jl")
 include("transfermatrix.jl")
 include("models/models.jl")
@@ -60,7 +60,7 @@ export Cell,
   InfMPS,
   InfiniteSum,
   InfiniteMPO,
-  InfiniteMPOMatrix,
+  InfiniteBlockMPO,
   InfiniteSumLocalOps,
   ITensorMap,
   ITensorNetwork,
diff --git a/src/infinitempo.jl b/src/infinitempo.jl
index 7a7fe69..7d0bcc4 100644
--- a/src/infinitempo.jl
+++ b/src/infinitempo.jl
@@ -60,11 +60,11 @@ function cat_to_itensor(Hm::Matrix{ITensor})
   return H, new_l, new_rs[1]
 end
 #
-#  Hm is the InfiniteMPOMatrix
+#  Hm is the InfiniteBlockMPO
 #  Hlrs is an array of {ITensor,Index,Index}s, one for each site in the unit cell.
 #  Hi is a CelledVector of ITensors.
 #
-function InfiniteMPO(Hm::InfiniteMPOMatrix)
+function InfiniteMPO(Hm::InfiniteBlockMPO)
   Hlrs = cat_to_itensor.(Hm) #return an array of {ITensor,Index,Index}
   #
   #  Unpack the array of tuples into three arrays.  And also get an array site indices.
@@ -115,8 +115,8 @@ function fuse_legs!(Hcl::InfiniteMPO, L, R)
       bottom => uniqueinds(bottom, dag(right_link));
       tags=[tags(right_link)],
     )
-    Hcl.data.data[j] = Hcl[j] * comb
-    Hcl.data.data[j + 1] = dag(comb) * Hcl[j + 1]
+    Hcl[j] = Hcl[j] * comb
+    Hcl[j + 1] = dag(comb) * Hcl[j + 1]
   end
   right_link = only(commoninds(Hcl[N], Hcl[N + 1]))
   right_link2 = translatecell(translator(Hcl), right_link, -1)
@@ -135,8 +135,8 @@ function fuse_legs!(Hcl::InfiniteMPO, L, R)
   )
   comb2 = translatecell(translator(Hcl), comb, -1)
   comb2_ind = translatecell(translator(Hcl), comb2, -1)
-  Hcl.data.data[N] = Hcl[N] * comb
-  Hcl.data.data[1] = dag(comb2) * Hcl[1]
+  Hcl[N] = Hcl[N] * comb
+  Hcl[1] = dag(comb2) * Hcl[1]
   L = L * comb2
   R = dag(comb) * R
   return L, R
diff --git a/src/infinitempomatrix.jl b/src/infinitempoblock.jl
similarity index 82%
rename from src/infinitempomatrix.jl
rename to src/infinitempoblock.jl
index a00678a..e8f0657 100644
--- a/src/infinitempomatrix.jl
+++ b/src/infinitempoblock.jl
@@ -1,16 +1,16 @@
 
-mutable struct InfiniteMPOMatrix <: AbstractInfiniteMPS
+mutable struct InfiniteBlockMPO <: AbstractInfiniteMPS
   data::CelledVector{Matrix{ITensor}}
   llim::Int #RealInfinity
   rlim::Int #RealInfinity
   reverse::Bool
 end
 
-translator(mpo::InfiniteMPOMatrix) = mpo.data.translator
-data(mpo::InfiniteMPOMatrix) = mpo.data
+translator(mpo::InfiniteBlockMPO) = mpo.data.translator
+data(mpo::InfiniteBlockMPO) = mpo.data
 
 # TODO better printing?
-function Base.show(io::IO, M::InfiniteMPOMatrix)
+function Base.show(io::IO, M::InfiniteBlockMPO)
   print(io, "$(typeof(M))")
   (length(M) > 0) && print(io, "\n")
   for i in eachindex(M)
@@ -19,8 +19,8 @@ function Base.show(io::IO, M::InfiniteMPOMatrix)
     else
       A = M[i]
       println(io, "Matrix tensor of size $(size(A))")
-      for k in 1:size(A)[1], l in 1:size(A)[2]
-        if !isassigned(A, k + (size(A)[1] - 1) * l)
+      for k in 1:size(A, 1), l in 1:size(A, 2)
+        if !isassigned(A, k + (size(A, 1) - 1) * l)
           println(io, "[($k, $l)] #undef")
         elseif isempty(A[k, l])
           println(io, "[($k, $l)] empty")
@@ -32,27 +32,27 @@ function Base.show(io::IO, M::InfiniteMPOMatrix)
   end
 end
 
-function getindex(ψ::InfiniteMPOMatrix, n::Integer)
+function getindex(ψ::InfiniteBlockMPO, n::Integer)
   return ψ.data[n]
 end
 
-function InfiniteMPOMatrix(arrMat::Vector{Matrix{ITensor}})
-  return InfiniteMPOMatrix(arrMat, 0, size(arrMat)[1], false)
+function InfiniteBlockMPO(arrMat::Vector{Matrix{ITensor}})
+  return InfiniteBlockMPO(arrMat, 0, size(arrMat, 1), false)
 end
 
-function InfiniteMPOMatrix(data::Vector{Matrix{ITensor}}, translator::Function)
-  return InfiniteMPOMatrix(CelledVector(data, translator), 0, size(data)[1], false)
+function InfiniteBlockMPO(data::Vector{Matrix{ITensor}}, translator::Function)
+  return InfiniteBlockMPO(CelledVector(data, translator), 0, size(data, 1), false)
 end
 
-function InfiniteMPOMatrix(data::CelledVector{Matrix{ITensor}}, m::Int64, n::Int64)
-  return InfiniteMPOMatrix(data, m, n, false)
+function InfiniteBlockMPO(data::CelledVector{Matrix{ITensor}}, m::Int64, n::Int64)
+  return InfiniteBlockMPO(data, m, n, false)
 end
 
-function InfiniteMPOMatrix(data::CelledVector{Matrix{ITensor}})
-  return InfiniteMPOMatrix(data, 0, size(data)[1], false)
+function InfiniteBlockMPO(data::CelledVector{Matrix{ITensor}})
+  return InfiniteBlockMPO(data, 0, size(data, 1), false)
 end
 
-function ITensors.siteinds(A::InfiniteMPOMatrix)
+function ITensors.siteinds(A::InfiniteBlockMPO)
   data = [dag(only(filterinds(uniqueinds(A[1][1, 1], A[2][1, 1]); plev=0)))]
   for x in 2:(nsites(A) - 1)
     append!(
@@ -71,7 +71,7 @@ function ITensors.siteinds(A::InfiniteMPOMatrix)
   return CelledVector(data, translator(A))
 end
 
-function ITensors.splitblocks(H::InfiniteMPOMatrix)
+function ITensors.splitblocks(H::InfiniteBlockMPO)
   H = copy(H)
   N = nsites(H)
   for j in 1:N
@@ -115,22 +115,10 @@ function find_all_links(Hm::Matrix{ITensor})
   return left_links, right_links
 end
 
-# function convert_itensor_to_itensormatrix(tensor; kwargs...)
-#   if order(tensor) == 3
-#     return convert_itensor_3vector(tensor; kwargs...)
-#   elseif order(tensor) == 4
-#     return convert_itensor_33matrix(tensor; kwargs...)
-#   else
-#     error(
-#       "Conversion of ITensor into matrix of ITensor not planned for this type of tensors"
-#     )
-#   end
-# end
-
 """
     local_mpo_block_projectors(is::Index; new_tags = tags(is))
 
-  Build the projectors on the three parts of the itensor used to split a MPO into an InfiniteMPOMatrix
+  Build the projectors on the three parts of the itensor used to split a MPO into an InfiniteBlockMPO
   More precisely, create projectors on the first dimension, the 2:end-1 and the last dimension of the index
   Input: is the Index to split
   Output: the triplet of projectors (first, middle, last)
diff --git a/src/models/models.jl b/src/models/models.jl
index fba98ab..36c84cd 100644
--- a/src/models/models.jl
+++ b/src/models/models.jl
@@ -87,14 +87,14 @@ function InfiniteMPO(model::Model, s::CelledVector; kwargs...)
 end
 
 function InfiniteMPO(model::Model, s::CelledVector, translator::Function; kwargs...)
-  return InfiniteMPO(InfiniteMPOMatrix(model, s, translator; kwargs...))
+  return InfiniteMPO(InfiniteBlockMPO(model, s, translator; kwargs...))
 end
 
-function InfiniteMPOMatrix(model::Model, s::CelledVector; kwargs...)
-  return InfiniteMPOMatrix(model, s, translator(s); kwargs...)
+function InfiniteBlockMPO(model::Model, s::CelledVector; kwargs...)
+  return InfiniteBlockMPO(model, s, translator(s); kwargs...)
 end
 
-function InfiniteMPOMatrix(model::Model, s::CelledVector, translator::Function; kwargs...)
+function InfiniteBlockMPO(model::Model, s::CelledVector, translator::Function; kwargs...)
   N = length(s)
   temp_H = InfiniteSum{MPO}(model, s; kwargs...)
   range_H = maximum(nrange(temp_H)) #Should be improved
@@ -168,7 +168,7 @@ function InfiniteMPOMatrix(model::Model, s::CelledVector, translator::Function;
     #mpos[j] += dense(Hmat) * setelt(ls[j-1] => total_dim) * setelt(ls[j] => total_dim)
   end
   #unify_indices and add virtual indices to the empty tensors
-  mpos = InfiniteMPOMatrix(mpos, translator)
+  mpos = InfiniteBlockMPO(mpos, translator)
   for x in 1:N
     sp = prime(s[x])
     sd = dag(s[x])
@@ -195,25 +195,25 @@ function InfiniteMPOMatrix(model::Model, s::CelledVector, translator::Function;
     for j in 2:(size(mpos[x], 1) - 1)
       for k in 2:(size(mpos[x], 2) - 1)
         if isempty(mpos[x][j, k])
-          mpos.data.data[x][j, k] = ITensor(left_inds[j - 1], sd, sp, new_right_inds[k - 1])
+          mpos[x][j, k] = ITensor(left_inds[j - 1], sd, sp, new_right_inds[k - 1])
         else
-          replaceinds!(mpos.data.data[x][j, k], right_inds[k - 1] => new_right_inds[k - 1])
+          replaceinds!(mpos[x][j, k], right_inds[k - 1] => new_right_inds[k - 1])
         end
       end
     end
     for j in [1, size(mpos[x], 1)]
       for k in 2:(size(mpos[x], 2) - 1)
         if isempty(mpos[x][j, k])
-          mpos.data.data[x][j, k] = ITensor(sd, sp, new_right_inds[k - 1])
+          mpos[x][j, k] = ITensor(sd, sp, new_right_inds[k - 1])
         else
-          replaceinds!(mpos.data.data[x][j, k], right_inds[k - 1] => new_right_inds[k - 1])
+          replaceinds!(mpos[x][j, k], right_inds[k - 1] => new_right_inds[k - 1])
         end
       end
     end
     for j in 2:(size(mpos[x], 1) - 1)
       for k in [1, size(mpos[x], 2)]
         if isempty(mpos[x][j, k])
-          mpos.data.data[x][j, k] = ITensor(left_inds[j - 1], sd, sp)
+          mpos[x][j, k] = ITensor(left_inds[j - 1], sd, sp)
         end
       end
     end
diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl
index 5d2eb20..f2c9a23 100644
--- a/src/subspace_expansion.jl
+++ b/src/subspace_expansion.jl
@@ -118,7 +118,7 @@ function generate_twobody_nullspace(
 end
 
 function generate_twobody_nullspace(
-  ψ::InfiniteCanonicalMPS, H::InfiniteMPOMatrix, b::Tuple{Int,Int}; atol=1e-2
+  ψ::InfiniteCanonicalMPS, H::InfiniteBlockMPO, b::Tuple{Int,Int}; atol=1e-2
 )
   n_1, n_2 = b
   L, _ = left_environment(H, ψ)
diff --git a/src/vumps_mpo.jl b/src/vumps_mpo.jl
index 5421bf4..4a7c561 100644
--- a/src/vumps_mpo.jl
+++ b/src/vumps_mpo.jl
@@ -1,7 +1,7 @@
 # Struct for use in linear system solver
 struct AOᴸ
   ψ::InfiniteCanonicalMPS
-  H::InfiniteMPOMatrix
+  H::InfiniteBlockMPO
   n::Int
 end
 
@@ -29,7 +29,7 @@ function (A::AOᴸ)(x)
 end
 
 function initialize_left_environment(
-  H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n::Int64; init_last=false
+  H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n::Int64; init_last=false
 )
   dₕ = size(H[n + 1], 1)
   sit = inds(H[n + 1][1, 1])
@@ -49,7 +49,7 @@ function initialize_left_environment(
 end
 
 function apply_local_left_transfer_matrix(
-  Lstart::Vector{ITensor}, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64
+  Lstart::Vector{ITensor}, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64
 )
   dₕ = length(Lstart)
   ψ′ = dag(ψ)'
@@ -67,7 +67,7 @@ end
 
 # apply the left transfer matrix at position n1 to the vector Lstart considering it at position m, adding to Ltarget
 function apply_local_left_transfer_matrix(
-  Lstart::ITensor, m::Int64, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64;
+  Lstart::ITensor, m::Int64, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64;
 )
   Ltarget = initialize_left_environment(H, ψ, n_1; init_last=false)
   for j in 1:m
@@ -78,7 +78,7 @@ end
 
 #apply the left transfer matrix n1:n1+nsites(ψ)-1
 function apply_left_transfer_matrix(
-  Lstart::ITensor, m::Int64, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64
+  Lstart::ITensor, m::Int64, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64
 )
   Ltarget = apply_local_left_transfer_matrix(Lstart, m, H, ψ, n_1)
   for j in 1:(nsites(ψ) - 1)
@@ -89,7 +89,7 @@ end
 
 # Also input C bond matrices to help compute the right fixed points
 # of ψ (R ≈ C * dag(C))
-function left_environment(H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS; tol=1e-10)
+function left_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e-10)
   N = nsites(H)
   @assert N == nsites(ψ)
 
@@ -151,7 +151,7 @@ end
 # Struct for use in linear system solver
 struct AOᴿ
   ψ::InfiniteCanonicalMPS
-  H::InfiniteMPOMatrix
+  H::InfiniteBlockMPO
   n::Int
 end
 
@@ -179,7 +179,7 @@ function (A::AOᴿ)(x)
 end
 
 function initialize_right_environment(
-  H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n::Int64; init_first=false
+  H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n::Int64; init_first=false
 )
   dₕ = size(H[n - 1], 2)
   sit = inds(H[n - 1][1, 1])
@@ -199,7 +199,7 @@ function initialize_right_environment(
 end
 
 function apply_local_right_transfer_matrix(
-  Lstart::Vector{ITensor}, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64
+  Lstart::Vector{ITensor}, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64
 )
   dₕ = length(Lstart)
   ψ′ = dag(ψ)'
@@ -218,7 +218,7 @@ end
 function apply_local_right_transfer_matrix(
   Lstart::ITensor,
   m::Int64,
-  H::InfiniteMPOMatrix,
+  H::InfiniteBlockMPO,
   ψ::InfiniteCanonicalMPS,
   n_1::Int64;
   reset=true,
@@ -236,7 +236,7 @@ end
 
 #apply the right transfer matrix n1:n1+nsites(ψ)-1
 function apply_right_transfer_matrix(
-  Lstart::ITensor, m::Int64, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64
+  Lstart::ITensor, m::Int64, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64
 )
   Ltarget = apply_local_right_transfer_matrix(Lstart, m, H, ψ, n_1)
   for j in 1:(nsites(ψ) - 1)
@@ -245,7 +245,7 @@ function apply_right_transfer_matrix(
   return Ltarget
 end
 
-function right_environment(H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS; tol=1e-10)
+function right_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e-10)
   N = nsites(H)
   @assert N == nsites(ψ)
 
@@ -303,7 +303,7 @@ function right_environment(H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS; tol=1
   return CelledVector(Rs), eᵣ[1]
 end
 
-function vumps(H::InfiniteMPOMatrix, ψ::InfiniteMPS; kwargs...)
+function vumps(H::InfiniteBlockMPO, ψ::InfiniteMPS; kwargs...)
   return vumps(H, orthogonalize(ψ, :); kwargs...)
 end
 
@@ -347,7 +347,7 @@ end
 
 function tdvp_iteration_sequential(
   solver::Function,
-  H::InfiniteMPOMatrix,
+  H::InfiniteBlockMPO,
   ψ::InfiniteCanonicalMPS;
   (ϵᴸ!)=fill(1e-15, nsites(ψ)),
   (ϵᴿ!)=fill(1e-15, nsites(ψ)),
@@ -409,7 +409,7 @@ end
 
 function tdvp_iteration_parallel(
   solver::Function,
-  H::InfiniteMPOMatrix,
+  H::InfiniteBlockMPO,
   ψ::InfiniteCanonicalMPS;
   (ϵᴸ!)=fill(1e-15, nsites(ψ)),
   (ϵᴿ!)=fill(1e-15, nsites(ψ)),
diff --git a/test/test_iMPOConversions.jl b/test/test_iMPOConversions.jl
index 1c0fefe..e0a19d9 100644
--- a/test/test_iMPOConversions.jl
+++ b/test/test_iMPOConversions.jl
@@ -2,7 +2,7 @@ using ITensors
 using ITensorInfiniteMPS
 using Test
 
-#With the new definition of InfiniteMPOMatrix, the MPO is better behaved, and hence we need to be a bit more careful
+#With the new definition of InfiniteBlockMPO, the MPO is better behaved, and hence we need to be a bit more careful
 function expect_over_unit_cell(ψ::InfiniteCanonicalMPS, h::InfiniteSum{MPO})
   s = siteinds(ψ)
   Ncell = nsites(h)
@@ -61,7 +61,7 @@ function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteMPO)
   return expect(ψ, terminate(h)) #defer to src/infinitecanonicalmps.jl
 end
 
-function generate_edges(h::InfiniteMPOMatrix)
+function generate_edges(h::InfiniteBlockMPO)
   Ncell = nsites(h)
   # left termination vector
   Ls = Matrix{ITensor}(undef, 1, size(h[1], 1))
@@ -82,7 +82,7 @@ function generate_edges(h::InfiniteMPOMatrix)
   return Ls, Rs
 end
 
-function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteMPOMatrix)
+function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteBlockMPO)
   Ncell = nsites(h)
   L, R = generate_edges(h)
   l = commoninds(ψ.AL[0], ψ.AL[1])
@@ -94,7 +94,7 @@ function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteMPOMatrix)
     temp = map(a -> contract(a, ψ.AR[j], dag(prime(ψ.AR[j]))), h[j])
     L = L * temp
   end
-  return (L * R)[1][1]#ITensorInfiniteMPS.scalar_product(L, R)[1]
+  return (L * R)[][]
 end
 
 #H = ΣⱼΣn (½ S⁺ⱼS⁻ⱼ₊n + ½ S⁻ⱼS⁺ⱼ₊n + SᶻⱼSᶻⱼ₊n)
@@ -139,7 +139,7 @@ function fermion_momentum_translator(i::Index, n::Integer; N=6)
   return new_i
 end
 
-@testset verbose = true "InfiniteMPOMatrix -> InfiniteMPO" begin
+@testset verbose = true "InfiniteBlockMPO -> InfiniteMPO" begin
   ferro(n) = "↑"
   antiferro(n) = isodd(n) ? "↑" : "↓"
 
@@ -159,7 +159,7 @@ end
     s = infsiteinds(site, Ncell; initstate, conserve_qns=qns)
     ψ = InfMPS(s, initstate)
     Hi = InfiniteMPO(model, s; model_kwargs...)
-    Hm = InfiniteMPOMatrix(model, s; model_kwargs...)
+    Hm = InfiniteBlockMPO(model, s; model_kwargs...)
     Hs = InfiniteSum{MPO}(model, s; model_kwargs...)
     Es = expect_over_unit_cell(ψ, Hs)
     Ei = expect(ψ, Hi)
diff --git a/test/test_vumpsmpo.jl b/test/test_vumpsmpo.jl
index c54c7cb..501790b 100644
--- a/test/test_vumpsmpo.jl
+++ b/test/test_vumpsmpo.jl
@@ -61,7 +61,7 @@ end
     s = infsiteinds("S=1/2", N; initstate, conserve_szparity=conserve_qns)
     ψ = InfMPS(s, initstate)
 
-    Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...)
+    Hmpo = InfiniteBlockMPO(model, s; model_kwargs...)
     # Alternate steps of running VUMPS and increasing the bond dimension
     ψ = tdvp(Hmpo, ψ; vumps_kwargs...)
     for _ in 1:outer_iters
@@ -133,7 +133,7 @@ end
     s = infsiteinds("S=1/2", nsites; conserve_szparity=conserve_qns, initstate)
     ψ = InfMPS(s, initstate)
 
-    Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...)
+    Hmpo = InfiniteBlockMPO(model, s; model_kwargs...)
     # Alternate steps of running VUMPS and increasing the bond dimension
     ψ = tdvp(Hmpo, ψ; vumps_kwargs...)
     for _ in 1:outer_iters
@@ -218,7 +218,7 @@ end
     )
     ψ = InfMPS(s, initstate)
 
-    Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...)
+    Hmpo = InfiniteBlockMPO(model, s; model_kwargs...)
     # Alternate steps of running VUMPS and increasing the bond dimension
     ψ = tdvp(Hmpo, ψ; vumps_kwargs...)
     for _ in 1:outer_iters
diff --git a/test/test_vumpsmpo_fqhe.jl b/test/test_vumpsmpo_fqhe.jl
index 98a156f..e940635 100644
--- a/test/test_vumpsmpo_fqhe.jl
+++ b/test/test_vumpsmpo_fqhe.jl
@@ -84,7 +84,7 @@ end
     )
     ψ = InfMPS(s, initstate)
 
-    Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...)
+    Hmpo = InfiniteBlockMPO(model, s; model_kwargs...)
     # Alternate steps of running VUMPS and increasing the bond dimension
     ψ = tdvp(Hmpo, ψ; vumps_kwargs...)
     for _ in 1:outer_iters