diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index d7fd0da07..4bab6bf82 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -173,6 +173,7 @@ include("algorithms/statmech/leading_boundary.jl")
 include("algorithms/statmech/vumps.jl")
 include("algorithms/statmech/vomps.jl")
 include("algorithms/statmech/gradient_grassmann.jl")
+include("algorithms/statmech/idmrg.jl")
 
 include("algorithms/fidelity_susceptibility.jl")
 
diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl
index 0ca3efda0..adb4aeec3 100644
--- a/src/algorithms/approximate/idmrg.jl
+++ b/src/algorithms/approximate/idmrg.jl
@@ -21,7 +21,7 @@ function approximate!(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili
             end
 
             # right to left sweep
-            for col in size(ψ, 2):-1:1
+            for col in reverse(1:size(ψ, 2))
                 for row in 1:size(ψ, 1)
                     ψ.AC[row + 1, col] = AC_projection(CartesianIndex(row, col),
                                                        ψ, toapprox, envs)
@@ -71,38 +71,42 @@ function approximate!(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili
             C_current = ψ.C[:, 0]
 
             # sweep from left to right
-            for col in 1:size(ψ, 2)
+            for site in 1:size(ψ, 2)
                 for row in 1:size(ψ, 1)
-                    AC2′ = AC2_projection(CartesianIndex(row, col), ψ, toapprox, envs)
+                    AC2′ = AC2_projection(CartesianIndex(row, site), ψ, toapprox, envs;
+                                          kind=:ACAR)
                     al, c, ar, = tsvd!(AC2′; trunc=alg.trscheme, alg=alg.alg_svd)
                     normalize!(c)
 
-                    ψ.AL[row + 1, col] = al
-                    ψ.C[row + 1, col] = complex(c)
-                    ψ.AR[row + 1, col + 1] = _transpose_front(ar)
-                    ψ.AC[row + 1, col + 1] = _transpose_front(c * ar)
+                    ψ.AL[row + 1, site] = al
+                    ψ.C[row + 1, site] = complex(c)
+                    ψ.AR[row + 1, site + 1] = _transpose_front(ar)
+                    ψ.AC[row + 1, site + 1] = _transpose_front(c * ar)
                 end
-                transfer_leftenv!(envs, ψ, toapprox, col + 1)
-                transfer_rightenv!(envs, ψ, toapprox, col)
+
+                transfer_leftenv!(envs, ψ, toapprox, site + 1)
+                transfer_rightenv!(envs, ψ, toapprox, site)
             end
+
             normalize!(envs, ψ, toapprox)
 
             # sweep from right to left
-            for col in (size(ψ, 2) - 1):-1:0
+            for site in reverse(0:(size(ψ, 2) - 1))
                 for row in 1:size(ψ, 1)
-                    # TODO: also write this as AC2_projection?
-                    AC2 = ϕ.AL[row, col] * _transpose_tail(ϕ.AC[row, col + 1])
-                    AC2′ = AC2_hamiltonian(CartesianIndex(row, col), ψ, O, ϕ, envs) * AC2
+                    AC2′ = AC2_projection(CartesianIndex(row, site), ψ, toapprox, envs;
+                                          kind=:ALAC)
                     al, c, ar, = tsvd!(AC2′; trunc=alg.trscheme, alg=alg.alg_svd)
                     normalize!(c)
 
-                    ψ.AL[row + 1, col] = al
-                    ψ.C[row + 1, col] = complex(c)
-                    ψ.AR[row + 1, col + 1] = _transpose_front(ar)
+                    ψ.AL[row + 1, site] = al
+                    ψ.C[row + 1, site] = complex(c)
+                    ψ.AR[row + 1, site + 1] = _transpose_front(ar)
                 end
-                transfer_leftenv!(envs, ψ, toapprox, col + 1)
-                transfer_rightenv!(envs, ψ, toapprox, col)
+
+                transfer_leftenv!(envs, ψ, toapprox, site + 1)
+                transfer_rightenv!(envs, ψ, toapprox, site)
             end
+
             normalize!(envs, ψ, toapprox)
 
             # update error
@@ -127,7 +131,7 @@ function approximate!(ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO,<:Multili
 
     # TODO: immediately compute in-place
     alg_gauge = updatetol(alg.alg_gauge, iter, ϵ)
-    ψ′ = MultilineMPS(map(x -> x, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
+    ψ′ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
     copy!(ψ, ψ′) # ensure output destination is unchanged
 
     recalculate!(envs, ψ, toapprox)
diff --git a/src/algorithms/derivatives/derivatives.jl b/src/algorithms/derivatives/derivatives.jl
index 294b029cf..869da571e 100644
--- a/src/algorithms/derivatives/derivatives.jl
+++ b/src/algorithms/derivatives/derivatives.jl
@@ -164,20 +164,21 @@ for kind in (:C, :AC, :AC2)
     projection = Symbol(kind, "_projection")
     hamiltonian = Symbol(kind, "_hamiltonian")
 
-    @eval function $projection(site, below, above::Tuple, envs)
-        return $projection(site, below, above..., envs)
+    @eval function $projection(site, below, above::Tuple, envs; kwargs...)
+        return $projection(site, below, above..., envs; kwargs...)
     end
-    @eval function $projection(site, below, above::AbstractMPS, envs)
-        return $projection(site, below, nothing, above, envs)
+    @eval function $projection(site, below, above::AbstractMPS, envs; kwargs...)
+        return $projection(site, below, nothing, above, envs; kwargs...)
     end
     @eval function $projection(site::CartesianIndex{2}, below::MultilineMPS, operator,
-                               above::MultilineMPS, envs)
+                               above::MultilineMPS, envs; kwargs...)
         row, col = Tuple(site)
-        return $projection(col, below[row + 1], operator[row], above[row], envs[row])
+        return $projection(col, below[row + 1], operator[row], above[row], envs[row];
+                           kwargs...)
     end
-    @eval function $projection(site, below, above::LazySum, envs)
+    @eval function $projection(site, below, above::LazySum, envs; kwargs...)
         return sum(zip(above.ops, envs.envs)) do x
-            return $projection(site, below, x...)
+            return $projection(site, below, x...; kwargs...)
         end
     end
 end
@@ -187,9 +188,8 @@ end
 function AC_projection(site, below, operator, above, envs)
     return AC_hamiltonian(site, below, operator, above, envs) * above.AC[site]
 end
-function AC2_projection(site::Int, below, operator, above, envs)
-    AC2 = above.AC[site] * _transpose_tail(above.AR[site + 1])
-    return AC2_hamiltonian(site, below, operator, above, envs) * AC2
+function AC2_projection(site::Int, below, operator, above, envs; kwargs...)
+    return AC2_hamiltonian(site, below, operator, above, envs) * AC2(above, site; kwargs...)
 end
 
 # Multiline
diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl
index afe81d055..0d0fd1843 100644
--- a/src/algorithms/groundstate/idmrg.jl
+++ b/src/algorithms/groundstate/idmrg.jl
@@ -40,7 +40,12 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG, envs=environments(ost
             for pos in 1:length(ψ)
                 h = AC_hamiltonian(pos, ψ, H, ψ, envs)
                 _, ψ.AC[pos] = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve)
-                ψ.AL[pos], ψ.C[pos] = leftorth!(ψ.AC[pos])
+                if pos == length(ψ)
+                    # AC needed in next sweep
+                    ψ.AL[pos], ψ.C[pos] = leftorth(ψ.AC[pos])
+                else
+                    ψ.AL[pos], ψ.C[pos] = leftorth!(ψ.AC[pos])
+                end
                 transfer_leftenv!(envs, ψ, H, ψ, pos + 1)
             end
 
@@ -124,7 +129,7 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(os
 
             # sweep from left to right
             for pos in 1:(length(ψ) - 1)
-                ac2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
+                ac2 = AC2(ψ, pos; kind=:ACAR)
                 h_ac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs)
                 _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve)
 
@@ -141,20 +146,22 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(os
             end
 
             # update the edge
-            @plansor ac2[-1 -2; -3 -4] := ψ.AC[end][-1 -2; 1] * inv(ψ.C[0])[1; 2] *
-                                          ψ.AL[1][2 -4; 3] * ψ.C[1][3; -3]
+            ψ.AL[end] = ψ.AC[end] / ψ.C[end]
+            ψ.AC[1] = _mul_tail(ψ.AL[1], ψ.C[1])
+            ac2 = AC2(ψ, 0; kind=:ALAC)
             h_ac2 = AC2_hamiltonian(0, ψ, H, ψ, envs)
             _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve)
 
             al, c, ar, = tsvd!(ac2′; trunc=alg.trscheme, alg=alg.alg_svd)
             normalize!(c)
 
-            ψ.AC[end] = al * c
             ψ.AL[end] = al
             ψ.C[end] = complex(c)
             ψ.AR[1] = _transpose_front(ar)
+
+            ψ.AC[end] = _mul_tail(al, c)
             ψ.AC[1] = _transpose_front(c * ar)
-            ψ.AL[1] = ψ.AC[1] * inv(ψ.C[1])
+            ψ.AL[1] = ψ.AC[1] / ψ.C[1]
 
             C_current = complex(c)
 
@@ -164,7 +171,7 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(os
 
             # sweep from right to left
             for pos in (length(ψ) - 1):-1:1
-                ac2 = ψ.AL[pos] * _transpose_tail(ψ.AC[pos + 1])
+                ac2 = AC2(ψ, pos; kind=:ALAC)
                 h_ac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs)
                 _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve)
 
@@ -172,7 +179,7 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(os
                 normalize!(c)
 
                 ψ.AL[pos] = al
-                ψ.AC[pos] = al * c
+                ψ.AC[pos] = _mul_tail(al, c)
                 ψ.C[pos] = complex(c)
                 ψ.AR[pos + 1] = _transpose_front(ar)
                 ψ.AC[pos + 1] = _transpose_front(c * ar)
@@ -182,17 +189,19 @@ function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs=environments(os
             end
 
             # update the edge
-            @plansor ac2[-1 -2; -3 -4] := ψ.C[end - 1][-1; 1] * ψ.AR[end][1 -2; 2] *
-                                          inv(ψ.C[end])[2; 3] * ψ.AC[1][3 -4; -3]
+            ψ.AC[end] = _mul_front(ψ.C[end - 1], ψ.AR[end])
+            ψ.AR[1] = _transpose_front(ψ.C[end] \ _transpose_tail(ψ.AC[1]))
+            ac2 = AC2(ψ, 0; kind=:ACAR)
             h_ac2 = AC2_hamiltonian(0, ψ, H, ψ, envs)
             _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve)
             al, c, ar, = tsvd!(ac2′; trunc=alg.trscheme, alg=alg.alg_svd)
             normalize!(c)
 
-            ψ.AR[end] = _transpose_front(inv(ψ.C[end - 1]) * _transpose_tail(al * c))
             ψ.AL[end] = al
             ψ.C[end] = complex(c)
             ψ.AR[1] = _transpose_front(ar)
+
+            ψ.AR[end] = _transpose_front(ψ.C[end - 1] \ _transpose_tail(al * c))
             ψ.AC[1] = _transpose_front(c * ar)
 
             transfer_leftenv!(envs, ψ, H, ψ, 1)
diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl
new file mode 100644
index 000000000..515cb805d
--- /dev/null
+++ b/src/algorithms/statmech/idmrg.jl
@@ -0,0 +1,194 @@
+function leading_boundary(ψ::MultilineMPS, operator, alg::IDMRG,
+                          envs=environments(ψ, operator))
+    log = IterLog("IDMRG")
+    ϵ::Float64 = 2 * alg.tol
+    local iter
+
+    LoggingExtras.withlevel(; alg.verbosity) do
+        @infov 2 loginit!(log, ϵ, expectation_value(ψ, operator, envs))
+        for outer iter in 1:(alg.maxiter)
+            alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ)
+            C_current = ψ.C[:, 0]
+
+            # left to right sweep
+            for col in 1:size(ψ, 2)
+                Hac = AC_hamiltonian(col, ψ, operator, ψ, envs)
+                _, ψ.AC[:, col] = fixedpoint(Hac, ψ.AC[:, col], :LM, alg_eigsolve)
+
+                for row in 1:size(ψ, 1)
+                    ac = ψ.AC[row, col]
+                    (col == size(ψ, 2)) && (ac = copy(ac)) # needed in next sweep
+                    ψ.AL[row, col], ψ.C[row, col] = leftorth!(ac)
+                end
+
+                transfer_leftenv!(envs, ψ, operator, ψ, col + 1)
+            end
+
+            # right to left sweep
+            for col in size(ψ, 2):-1:1
+                Hac = AC_hamiltonian(col, ψ, operator, ψ, envs)
+                _, ψ.AC[:, col] = fixedpoint(Hac, ψ.AC[:, col], :LM, alg_eigsolve)
+
+                for row in 1:size(ψ, 1)
+                    ψ.C[row, col - 1], temp = rightorth!(_transpose_tail(ψ.AC[row, col]))
+                    ψ.AR[row, col] = _transpose_front(temp)
+                end
+
+                transfer_rightenv!(envs, ψ, operator, ψ, col - 1)
+            end
+
+            normalize!(envs, ψ, operator, ψ)
+
+            ϵ = norm(C_current - ψ.C[:, 0])
+
+            if ϵ < alg.tol
+                @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, operator, envs))
+                break
+            end
+            if iter == alg.maxiter
+                @warnv 1 logcancel!(log, iter, ϵ, expectation_value(ψ, operator, envs))
+            else
+                @infov 3 logiter!(log, iter, ϵ, expectation_value(ψ, operator, envs))
+            end
+        end
+    end
+
+    alg_gauge = updatetol(alg.alg_gauge, iter, ϵ)
+    ψ = MultilineMPS(map(x -> x, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
+
+    recalculate!(envs, ψ, operator, ψ)
+    return ψ, envs, ϵ
+end
+
+function leading_boundary(ψ::MultilineMPS, operator, alg::IDMRG2,
+                          envs=environments(ψ, operator))
+    size(ψ, 2) < 2 && throw(ArgumentError("unit cell should be >= 2"))
+    ϵ::Float64 = 2 * alg.tol
+    log = IterLog("IDMRG2")
+    local iter
+
+    LoggingExtras.withlevel(; alg.verbosity) do
+        @infov 2 loginit!(log, ϵ)
+        for outer iter in 1:(alg.maxiter)
+            alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ)
+            C_current = ψ.C[:, 0]
+
+            # sweep from left to right
+            for site in 1:(size(ψ, 2) - 1)
+                ac2 = AC2(ψ, site; kind=:ACAR)
+                h = AC2_hamiltonian(site, ψ, operator, ψ, envs)
+                _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
+
+                for row in 1:size(ψ, 1)
+                    al, c, ar, = tsvd!(ac2′[row]; trunc=alg.trscheme, alg=alg.alg_svd)
+                    normalize!(c)
+
+                    ψ.AL[row + 1, site] = al
+                    ψ.C[row + 1, site] = complex(c)
+                    ψ.AR[row + 1, site + 1] = _transpose_front(ar)
+                    ψ.AC[row + 1, site + 1] = _transpose_front(c * ar)
+                end
+
+                transfer_leftenv!(envs, ψ, operator, ψ, site + 1)
+                transfer_rightenv!(envs, ψ, operator, ψ, site)
+            end
+
+            normalize!(envs, ψ, operator, ψ)
+
+            # update the edge
+            site = size(ψ, 2)
+            ψ.AL[:, end] .= ψ.AC[:, end] ./ ψ.C[:, end]
+            ψ.AC[:, 1] .= _mul_tail.(ψ.AL[:, 1], ψ.C[:, 1])
+            ac2 = AC2(ψ, site; kind=:ALAC)
+            h = AC2_hamiltonian(site, ψ, operator, ψ, envs)
+            _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
+
+            for row in 1:size(ψ, 1)
+                al, c, ar, = tsvd!(ac2′[row]; trunc=alg.trscheme, alg=alg.alg_svd)
+                normalize!(c)
+
+                ψ.AL[row + 1, site] = al
+                ψ.C[row + 1, site] = complex(c)
+                ψ.AR[row + 1, site + 1] = _transpose_front(ar)
+
+                ψ.AC[row + 1, site] = _mul_tail(al, c)
+                ψ.AC[row + 1, 1] = _transpose_front(c * ar)
+                ψ.AL[row + 1, 1] = ψ.AC[row + 1, 1] / ψ.C[row + 1, 1]
+            end
+
+            # TODO: decide if we should compare at the half-sweep level?
+            # C_current = ψ.C[:, site]
+
+            transfer_leftenv!(envs, ψ, operator, ψ, 1)
+            transfer_rightenv!(envs, ψ, operator, ψ, 0)
+
+            # sweep from right to left
+            for site in reverse(1:(size(ψ, 2) - 1))
+                ac2 = AC2(ψ, site; kind=:ALAC)
+                h = AC2_hamiltonian(site, ψ, operator, ψ, envs)
+                _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
+
+                for row in 1:size(ψ, 1)
+                    al, c, ar, = tsvd!(ac2′[row]; trunc=alg.trscheme, alg=alg.alg_svd)
+                    normalize!(c)
+
+                    ψ.AL[row + 1, site] = al
+                    ψ.C[row + 1, site] = complex(c)
+                    ψ.AR[row + 1, site + 1] = _transpose_front(ar)
+                end
+
+                transfer_leftenv!(envs, ψ, operator, ψ, site + 1)
+                transfer_rightenv!(envs, ψ, operator, ψ, site)
+            end
+
+            normalize!(envs, ψ, operator, ψ)
+
+            # update the edge
+            ψ.AC[:, end] .= _mul_front.(ψ.C[:, end - 1], ψ.AR[:, end])
+            ψ.AR[:, 1] .= _transpose_front.(ψ.C[:, end] .\ _transpose_tail.(ψ.AC[:, 1]))
+            ac2 = AC2(ψ, 0; kind=:ACAR)
+            h = AC2_hamiltonian(0, ψ, operator, ψ, envs)
+            _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
+
+            for row in 1:size(ψ, 1)
+                al, c, ar, = tsvd!(ac2′[row]; trunc=alg.trscheme, alg=alg.alg_svd)
+                normalize!(c)
+
+                ψ.AL[row + 1, end] = al
+                ψ.C[row + 1, end] = complex(c)
+                ψ.AR[row + 1, 1] = _transpose_front(ar)
+
+                ψ.AR[row + 1, end] = _transpose_front(ψ.C[row + 1, end - 1] \
+                                                      _transpose_tail(al * c))
+                ψ.AC[row + 1, 1] = _transpose_front(c * ar)
+            end
+
+            transfer_leftenv!(envs, ψ, operator, ψ, 1)
+            transfer_rightenv!(envs, ψ, operator, ψ, 0)
+
+            # update error
+            ϵ = sum(zip(C_current, ψ.C[:, 0])) do (c1, c2)
+                smallest = infimum(_firstspace(c1), _firstspace(c2))
+                e1 = isometry(_firstspace(c1), smallest)
+                e2 = isometry(_firstspace(c2), smallest)
+                return norm(e2' * c2 * e2 - e1' * c1 * e1)
+            end
+
+            if ϵ < alg.tol
+                @infov 2 logfinish!(log, iter, ϵ)
+                break
+            end
+            if iter == alg.maxiter
+                @warnv 1 logcancel!(log, iter, ϵ)
+            else
+                @infov 3 logiter!(log, iter, ϵ)
+            end
+        end
+    end
+
+    alg_gauge = updatetol(alg.alg_gauge, iter, ϵ)
+    ψ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
+
+    recalculate!(envs, ψ, operator, ψ)
+    return ψ, envs, ϵ
+end
diff --git a/src/environments/infinite_envs.jl b/src/environments/infinite_envs.jl
index 2e56a51d1..2f301a2b1 100644
--- a/src/environments/infinite_envs.jl
+++ b/src/environments/infinite_envs.jl
@@ -99,9 +99,15 @@ function compute_rightenvs!(envs::InfiniteEnvironments, below::InfiniteMPS,
     return λ, envs
 end
 
+# normalization convention of the environments:
+# - normalize the right environment to have norm 1
+# - normalize the left environment to have overlap 1
+# this avoids catastrophic blow-up of norms, while keeping the total normalized
+# and does not lead to issues for negative overlaps and real entries.
 function TensorKit.normalize!(envs::InfiniteEnvironments, below::InfiniteMPS,
                               operator::InfiniteMPO, above::InfiniteMPS)
     for i in 1:length(operator)
+        normalize!(envs.GRs[i])
         Hc = C_hamiltonian(i, below, operator, above, envs)
         λ = dot(below.C[i], Hc * above.C[i])
         scale!(envs.GLs[i + 1], inv(λ))
diff --git a/src/states/abstractmps.jl b/src/states/abstractmps.jl
index 851c10ce4..88d57404e 100644
--- a/src/states/abstractmps.jl
+++ b/src/states/abstractmps.jl
@@ -126,6 +126,17 @@ function makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg=QRpos())
     return A
 end
 
+# Tensor accessors
+# ----------------
+@doc """
+    AC2(ψ::AbstractMPS, i; kind=:ACAR)
+
+Obtain the two-site (center) gauge tensor at site `i` of the MPS `ψ`.
+If this hasn't been computed before, this can be computed as:
+- `kind=:ACAR` : AC[i] * AR[i+1]
+- `kind=:ALAC` : AL[i] * AC[i+1]
+""" AC2
+
 #===========================================================================================
 MPS types
 ===========================================================================================#
diff --git a/src/states/finitemps.jl b/src/states/finitemps.jl
index 3d306d3f9..b0eef3bcb 100644
--- a/src/states/finitemps.jl
+++ b/src/states/finitemps.jl
@@ -318,6 +318,9 @@ Base.@propagate_inbounds function Base.getindex(ψ::FiniteMPS, i::Int)
     end
 end
 
+# TODO: check where gauge center is to determine efficient kind
+AC2(psi::FiniteMPS, site::Int) = psi.AC[site] * _transpose_tail(psi.AR[site + 1])
+
 _complex_if_not_missing(x) = ismissing(x) ? x : complex(x)
 function Base.complex(mps::FiniteMPS)
     scalartype(mps) <: Complex && return mps
diff --git a/src/states/infinitemps.jl b/src/states/infinitemps.jl
index cf4c677c4..af3820384 100644
--- a/src/states/infinitemps.jl
+++ b/src/states/infinitemps.jl
@@ -212,6 +212,16 @@ end
 Utility
 ===========================================================================================#
 
+function AC2(ψ::InfiniteMPS, i::Integer; kind=:ACAR)
+    if kind == :ACAR
+        return ψ.AC[i] * _transpose_tail(ψ.AR[i + 1])
+    elseif kind == :ALAC
+        return ψ.AL[i] * _transpose_tail(ψ.AC[i + 1])
+    else
+        throw(ArgumentError("Invalid kind: $kind"))
+    end
+end
+
 Base.size(ψ::InfiniteMPS, args...) = size(ψ.AL, args...)
 Base.length(ψ::InfiniteMPS) = length(ψ.AL)
 Base.eltype(ψ::InfiniteMPS) = eltype(typeof(ψ))
diff --git a/src/states/multilinemps.jl b/src/states/multilinemps.jl
index 9e4e4b87d..6211b144d 100644
--- a/src/states/multilinemps.jl
+++ b/src/states/multilinemps.jl
@@ -57,6 +57,15 @@ function Base.getproperty(psi::MultilineMPS, prop::Symbol)
     end
 end
 
+function AC2(psi::MultilineMPS, site::CartesianIndex{2}; kwargs...)
+    return AC2(psi[site[1]], site[2]; kwargs...)
+end
+function AC2(psi::MultilineMPS, site::Int; kwargs...)
+    return map(1:size(psi, 1)) do row
+        return AC2(psi, CartesianIndex(row, site); kwargs...)
+    end
+end
+
 function Base.propertynames(::MultilineMPS)
     return (:AL, :AR, :AC, :C)
 end
diff --git a/test/algorithms.jl b/test/algorithms.jl
index caa6bdda0..4a10af9d3 100644
--- a/test/algorithms.jl
+++ b/test/algorithms.jl
@@ -430,18 +430,28 @@ end
 @testset "leading_boundary" verbose = true begin
     tol = 1e-4
     verbosity = verbosity_conv
+    D = 10
+    D1 = 13
     algs = [VUMPS(; tol, verbosity), VOMPS(; tol, verbosity),
-            GradientGrassmann(; tol, verbosity)]
+            GradientGrassmann(; tol, verbosity), IDMRG(; tol, verbosity),
+            IDMRG2(; tol, verbosity, trscheme=truncdim(D1))]
     mpo = force_planar(classical_ising())
 
-    ψ₀ = InfiniteMPS([ℙ^2], [ℙ^10])
+    ψ₀ = InfiniteMPS([ℙ^2], [ℙ^D])
     @testset "Infinite $i" for (i, alg) in enumerate(algs)
-        ψ, envs = leading_boundary(ψ₀, mpo, alg)
-        ψ, envs = changebonds(ψ, mpo, OptimalExpand(; trscheme=truncdim(3)), envs)
-        ψ, envs = leading_boundary(ψ, mpo, alg)
-
-        @test dim(space(ψ.AL[1, 1], 1)) == dim(space(ψ₀.AL[1, 1], 1)) + 3
-        @test expectation_value(ψ, mpo, envs) ≈ 2.5337 atol = 1e-3
+        if alg isa IDMRG2
+            ψ2 = repeat(ψ₀, 2)
+            mpo2 = repeat(mpo, 2)
+            ψ, envs = leading_boundary(ψ2, mpo2, alg)
+            @test dim(space(ψ.AL[1, 1], 1)) == dim(space(ψ₀.AL[1, 1], 1)) + (D1 - D)
+            @test expectation_value(ψ, mpo2, envs) ≈ 2.5337^2 atol = 1e-3
+        else
+            ψ, envs = leading_boundary(ψ₀, mpo, alg)
+            ψ, envs = changebonds(ψ, mpo, OptimalExpand(; trscheme=truncdim(D1 - D)), envs)
+            ψ, envs = leading_boundary(ψ, mpo, alg)
+            @test dim(space(ψ.AL[1, 1], 1)) == dim(space(ψ₀.AL[1, 1], 1)) + (D1 - D)
+            @test expectation_value(ψ, mpo, envs) ≈ 2.5337 atol = 1e-3
+        end
     end
 end