diff --git a/.gitignore b/.gitignore
index d2e3e08d..1fdef104 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,3 +3,4 @@ Manifest.toml
 .vscode
 docs/build/
 .DS_Store
+CLAUDE.md
diff --git a/README.md b/README.md
index 47891018..9c88923a 100644
--- a/README.md
+++ b/README.md
@@ -81,7 +81,7 @@ scheme = BTRG(T) # Bond-weighted TRG (excellent choice)
 data = run!(scheme, truncrank(16), maxiter(25)) # max bond-dimension of 16, for 25 iterations
 ```
 
-`data` now contains 26 norms of the tensor, 1 for every time the tensor was normalized. (By default there is a normalization step before the first coarse-graining step wich can be turned off by changing the kwarg `run!(...; finalize_beginning=false)`)
+`data` now contains 25 norms of the tensor, 1 for every time the tensor was normalized.
 
 Using these norms you could, for example, calculate the free energy of the critical classical Ising model:
 ```Julia
diff --git a/docs/src/finalizers.md b/docs/src/finalizers.md
index a4fa12d3..37a012fc 100644
--- a/docs/src/finalizers.md
+++ b/docs/src/finalizers.md
@@ -1,5 +1,5 @@
 # Finalizers
-At the end of every TNR step (and before the first step if `finalize_beginning=true` is chose in the `run!` function, which is default behaviour), the state of the scheme is finalized.
+At the end of every TNR step, the state of the scheme is finalized.
 
 By default this finalization process is as follow:
 
diff --git a/docs/src/index.md b/docs/src/index.md
index 0cd019e5..3e72fb44 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -58,7 +58,7 @@ T = classical_ising(ising_βc) # partition function of classical Ising model at
 scheme = BTRG(T) # Bond-weighted TRG (excellent choice)
 data = run!(scheme, truncrank(16), maxiter(25)) # max bond-dimension of 16, for 25 iterations
 ```
-`data` now contains 26 norms of the tensor, 1 for every time the tensor was normalized. (By default there is a normalization step before the first coarse-graining step wich can be turned off by changing the kwarg `run!(...; finalize_beginning=false)`)
+`data` now contains 25 norms of the tensor, 1 for every time the tensor was normalized.
 
 Using these norms you could, for example, calculate the free energy of the critical classical Ising model:
 ```Julia
diff --git a/examples/example.jl b/examples/example.jl
index 17c2bd11..4cf699aa 100644
--- a/examples/example.jl
+++ b/examples/example.jl
@@ -9,12 +9,8 @@ stopping_criterion = convcrit(1.0e-16, trg_f) & maxiter(20)
 # choose a TensorKit truncation scheme
 trunc = truncrank(16) & trunctol(atol = 1.0e-40)
 
-# initialize the TRG scheme
-scheme = TRG(classical_ising(1.0))
 
-# run the TRG scheme (and normalize and store the norm in the beginning (finalize_beginning=true))
-data = run!(scheme, trunc, stopping_criterion; finalize_beginning = true)
-# or: data = run!(scheme, truncrank(16)), this will default to maxiter(100)
+# ---- old interface (other schemes) ----
 
 # initialize the BTRG scheme
 scheme = BTRG(classical_ising(1.0), -0.5)
@@ -27,3 +23,30 @@ scheme = HOTRG(classical_ising(1.0))
 
 # run the HOTRG scheme
 data = run!(scheme, trunc, stopping_criterion)
+
+# ---- iterable `Renormalizer` interface ----
+
+# create a pure algorithm config (kwargs with sensible defaults)
+alg = TRG(; trunc = trunc, stop = stopping_criterion)
+renorm = Renormalizer(alg, classical_ising(1.0))
+
+# run all steps with logging
+state, data = run!(renorm; verbosity = 1)
+f = free_energy(data, 1.0)
+
+# or, step through manually
+renorm = Renormalizer(alg, classical_ising(1.0))
+for (state, data) in renorm
+    τ0, _ = extract_tau_and_c(state.T; fast = true)
+    # each iteration yields the state after one RG step
+    # renorm.step tracks how many steps have been completed
+end
+# renorm.step is the number of completed steps
+
+# or, advance one step at a time
+renorm = Renormalizer(alg, classical_ising(1.0))
+rgstep!(renorm)  # throws if stop criterion already met
+rgstep!(renorm)
+
+# access the current tensor directly from the state
+T_final = renorm.state.T
diff --git a/src/TNRKit.jl b/src/TNRKit.jl
index 78f47824..f7945139 100644
--- a/src/TNRKit.jl
+++ b/src/TNRKit.jl
@@ -24,6 +24,7 @@ export trivial_convcrit
 
 # schemes
 include("schemes/tnrscheme.jl")
+include("schemes/renormalizer.jl")
 include("schemes/trg.jl")
 include("schemes/btrg.jl")
 include("schemes/hotrg.jl")
@@ -61,9 +62,9 @@ include("schemes/looptnr.jl")
 include("schemes/symmetric_looptnr.jl")
 export classical_ising_inv # Ising model with all legs in codomain
 
-export TNRScheme
+export Renormalizer, TNRScheme, TNRAlgorithm
 
-export TRG
+export TRG, OneSiteState
 export BTRG
 export HOTRG
 export HOTRG_3D
@@ -91,7 +92,7 @@ export CorrelationHOTRG
 export LoopTNR, LoopParameters
 export SLoopTNR
 
-export run!
+export run!, rgstep!
 
 # models
 include("models/ising.jl")
diff --git a/src/schemes/atrg.jl b/src/schemes/atrg.jl
index a07b7c46..87df95a4 100644
--- a/src/schemes/atrg.jl
+++ b/src/schemes/atrg.jl
@@ -7,7 +7,7 @@ Anisotropic Tensor Renormalization Group
     $(FUNCTIONNAME)(T)
 
 # Running the algorithm
-    run!(::ATRG, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=default_Finalizer, finalize_beginning=true, verbosity=1])
+    run!(::ATRG, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=default_Finalizer, verbosity=1])
 
 Each step rescales the lattice by a (linear) factor of √2
 
diff --git a/src/schemes/atrg3d.jl b/src/schemes/atrg3d.jl
index bbf7843f..20f8d7bf 100644
--- a/src/schemes/atrg3d.jl
+++ b/src/schemes/atrg3d.jl
@@ -7,7 +7,7 @@ $(TYPEDEF)
     $(FUNCTIONNAME)(T)
 
 # Running the algorithm
-    run!(::ATRG_3D, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=defualt_Finalizer, finalize_beginning=true,verbosity=1])
+    run!(::ATRG_3D, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=defualt_Finalizer, verbosity=1])
 
 Each step rescales the lattice by a (linear) factor of 2
 
diff --git a/src/schemes/btrg.jl b/src/schemes/btrg.jl
index c2d4febc..f64a1366 100644
--- a/src/schemes/btrg.jl
+++ b/src/schemes/btrg.jl
@@ -7,7 +7,7 @@ Bond-weighted Tensor Renormalization Group
     $(FUNCTIONNAME)(T [, k=-1/2])
 
 # Running the algorithm
-    run!(::BTRG, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=default_Finalizer, finalize_beginning=true, verbosity=1])
+    run!(::BTRG, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=default_Finalizer, verbosity=1])
 
 Each step rescales the lattice by a (linear) factor of √2
 
diff --git a/src/schemes/ctm/c4vctm.jl b/src/schemes/ctm/c4vctm.jl
index e0af75bb..2cd38041 100644
--- a/src/schemes/ctm/c4vctm.jl
+++ b/src/schemes/ctm/c4vctm.jl
@@ -12,7 +12,7 @@ or with a (0,4) tensor (North, East, South, West) (unflipped arrow convention).
 The keyword argument symmetrize makes the tensor C4v symmetric when set to true. If symmetrize = false, it checks the symmetry explicitly.
 
 # Running the algorithm
-    run!(::c4vCTM, trunc::TruncationStrategy, stop::Stopcrit[, finalize_beginning=true, verbosity=1])
+    run!(::c4vCTM, trunc::TruncationStrategy, stop::Stopcrit[, verbosity=1])
 
 !!! info "verbosity levels"
     - 0: No output
diff --git a/src/schemes/ctm/honeycomb.jl b/src/schemes/ctm/honeycomb.jl
index d9124547..eb580362 100644
--- a/src/schemes/ctm/honeycomb.jl
+++ b/src/schemes/ctm/honeycomb.jl
@@ -25,7 +25,7 @@ or with a (0,3) tensor (120°, 0°, 240°) where all arrows point inward (unflip
 The keyword argument symmetrize makes the tensor C6v symmetric when set to true. If symmetrize = false, it checks the symmetry explicitly.
 
 # Running the algorithm
-    run!(::CTM, trunc::TruncationStrategy, stop::Stopcrit[, finalize_beginning=true, verbosity=1])
+    run!(::CTM, trunc::TruncationStrategy, stop::Stopcrit[, verbosity=1])
 
 !!! info "verbosity levels"
     - 0: No output
@@ -131,7 +131,7 @@ In the flipped arrow convention, the arrows point from (120°) to (240°, 0°).
 or with a (0,3) tensor (120°, 0°, 240°) where all arrows point inward (unflipped arrow convention).
 
 # Running the algorithm
-    run!(::CTM, trunc::TruncationStrategy, stop::Stopcrit[, finalize_beginning=true, verbosity=1])
+    run!(::CTM, trunc::TruncationStrategy, stop::Stopcrit[, verbosity=1])
 
 !!! info "verbosity levels"
     - 0: No output
diff --git a/src/schemes/ctm/triangular.jl b/src/schemes/ctm/triangular.jl
index 46e2dbfe..9c337c10 100644
--- a/src/schemes/ctm/triangular.jl
+++ b/src/schemes/ctm/triangular.jl
@@ -25,7 +25,7 @@ or with a (0,6) tensor (120°, 60°, 0°, 300°, 240°, 180°) where all arrows
 The keyword argument symmetrize makes the tensor C6v symmetric when set to true. If symmetrize = false, it checks the symmetry explicitly.
 
 # Running the algorithm
-    run!(::CTM, trunc::TruncationStrategy, stop::Stopcrit[, finalize_beginning=true, verbosity=1])
+    run!(::CTM, trunc::TruncationStrategy, stop::Stopcrit[, verbosity=1])
 
 !!! info "verbosity levels"
     - 0: No output
@@ -100,7 +100,7 @@ or with a (0,6) tensor (120°, 60°, 0°, 300°, 240°, 180°) where all arrows
 The keyword argument symmetrize makes the tensor C6v symmetric when set to true. If symmetrize = false, it checks the symmetry explicitly.
 
 # Running the algorithm
-    run!(::c6vCTM, trunc::MatrixAlgebraKit.TruncationStrategy, stop::Stopcrit[, finalize_beginning=true, projectors=:twothirds, conditioning=true, verbosity=1])
+    run!(::c6vCTM, trunc::MatrixAlgebraKit.TruncationStrategy, stop::Stopcrit[, projectors=:twothirds, conditioning=true, verbosity=1])
 
 `projectors` can either be :twothirds or :full, determining the type of projectors used in the renormalization step. This is based on https://arxiv.org/abs/2510.04907v1.
 `conditioning` determines whether to condition the second projector construction. This is based on https://doi.org/10.1103/PhysRevB.98.235148.
diff --git a/src/schemes/hotrg.jl b/src/schemes/hotrg.jl
index c6d787db..06b772ed 100644
--- a/src/schemes/hotrg.jl
+++ b/src/schemes/hotrg.jl
@@ -7,7 +7,7 @@ Higher-Order Tensor Renormalization Group
     $(FUNCTIONNAME)(T)
 
 # Running the algorithm
-    run!(::HOTRG, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=default_Finalizer, finalize_beginning=true, verbosity=1])
+    run!(::HOTRG, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=default_Finalizer, verbosity=1])
 
 Each step rescales the lattice by a (linear) factor of 2
 
diff --git a/src/schemes/hotrg3d.jl b/src/schemes/hotrg3d.jl
index 0be920a1..d2eabf94 100644
--- a/src/schemes/hotrg3d.jl
+++ b/src/schemes/hotrg3d.jl
@@ -7,7 +7,7 @@ $(TYPEDEF)
     $(FUNCTIONNAME)(T)
 
 # Running the algorithm
-    run!(::HOTRG_3D, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=default_Finalizer, finalize_beginning=true, verbosity=1])
+    run!(::HOTRG_3D, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=default_Finalizer, verbosity=1])
 
 Each step rescales the lattice by a (linear) factor of 2
 
diff --git a/src/schemes/impurityhotrg.jl b/src/schemes/impurityhotrg.jl
index 4d3b32ba..f38e8e4c 100644
--- a/src/schemes/impurityhotrg.jl
+++ b/src/schemes/impurityhotrg.jl
@@ -7,7 +7,7 @@ Single impurity method for Higher-Order Tensor Renormalization Group (for 2nd or
     $(FUNCTIONNAME)(T, T_imp_order1_1, T_imp_order1_2, T_imp_order2)
 
 # Running the algorithm
-    run!(::ImpurityHOTRG, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=ImpurityHOTRG_Finalizer, finalize_beginning=true, verbosity=1])
+    run!(::ImpurityHOTRG, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=ImpurityHOTRG_Finalizer, verbosity=1])
 
 Each step rescales the lattice by a (linear) factor of 2
 
diff --git a/src/schemes/impuritytrg.jl b/src/schemes/impuritytrg.jl
index 12235d11..4974b65b 100644
--- a/src/schemes/impuritytrg.jl
+++ b/src/schemes/impuritytrg.jl
@@ -7,7 +7,7 @@ Impurity method for Tensor Renormalization Group
     $(FUNCTIONNAME)(T, T_imp1, T_imp2, T_imp3, T_imp4)
 
 # Running the algorithm
-    run!(::ImpurityTRG, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=ImpurityTRG_Finalizer, finalize_beginning=true, verbosity=1])
+    run!(::ImpurityTRG, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=ImpurityTRG_Finalizer, verbosity=1])
 
 Each step rescales the lattice by a (linear) factor of √2
 
diff --git a/src/schemes/looptnr.jl b/src/schemes/looptnr.jl
index 8905cd94..6bf9b054 100644
--- a/src/schemes/looptnr.jl
+++ b/src/schemes/looptnr.jl
@@ -10,7 +10,7 @@ Loop Optimization for Tensor Network Renormalization
 
 # Running the algorithm
     run!(::LoopTNR, trunc::TruncationStrategy, criterion::stopcrit, [parameters::LoopParameters], [finalizer::Finalizer];
-        [entanglement_criterion::stopcrit, finalize_beginning=true, verbosity=1])
+        [entanglement_criterion::stopcrit, verbosity=1])
 
 # LoopParameters
 See also: [`LoopParameters`](@ref)
@@ -490,17 +490,12 @@ function run!(
         criterion::stopcrit, loop_condition::LoopParameters,
         finalizer::Finalizer{E};
         entanglement_criterion = default_entanglement_criterion,
-        finalize_beginning = true,
         verbosity = 1
     ) where {E}
     data = Vector{E}()
 
     LoggingExtras.withlevel(; verbosity) do
         @infov 1 "Starting simulation\n $(scheme)\n"
-        if finalize_beginning
-            push!(data, finalizer.f!(scheme))
-        end
-
         steps = 0
         crit = true
 
diff --git a/src/schemes/renormalizer.jl b/src/schemes/renormalizer.jl
new file mode 100644
index 00000000..8755b892
--- /dev/null
+++ b/src/schemes/renormalizer.jl
@@ -0,0 +1,77 @@
+abstract type TNRAlgorithm end
+
+"""
+$(TYPEDEF)
+
+Iterator that performs RG coarse-graining steps on network tensors.
+
+# Fields
+$(TYPEDFIELDS)
+
+# Iterator
+In a for-loop, each iteration yields `(state, data)` after each RG step.
+Stops when the algorithm's stopping criterion is met.
+"""
+mutable struct Renormalizer{A <: TNRAlgorithm, S, D}
+    "Algorithm configuration (truncation, maxiter, etc.)"
+    alg::A
+    "Algorithm-specific state holding all involved tensors"
+    state::S
+    "Accumulated per-step data"
+    data::Vector{D}
+    "Number of RG steps performed so far"
+    step::Int
+end
+
+Renormalizer(alg::A, state::S, data::Vector{D}) where {A, S, D} = Renormalizer{A, S, D}(alg, state, data, 0)
+Renormalizer(alg::A, state::S, data::Vector{D}, step::Int) where {A, S, D} = Renormalizer{A, S, D}(alg, state, data, step)
+
+function Base.iterate(r::Renormalizer)
+    return iterate(r, 0)
+end
+
+function Base.iterate(r::Renormalizer, n::Int)
+    r.alg.stop(n, r.data) || return nothing
+    step!(r.state, r.alg)
+    val = finalize!(r.state)
+    push!(r.data, val)
+    r.step = n + 1
+    return ((r.state, r.data), n + 1)
+end
+
+function Base.show(io::IO, r::Renormalizer)
+    println(io, "Renormalizer")
+    println(io, "  * algorithm: $(nameof(typeof(r.alg)))")
+    println(io, "  * state: $(nameof(typeof(r.state)))")
+    println(io, "  * step: $(r.step)")
+    println(io, "  * data: $(length(r.data)) entries")
+    return nothing
+end
+
+"""
+    rgstep!(r::Renormalizer)
+
+Perform one RG coarse-graining step. Wraps `Base.iterate`.
+Throws an error if the stopping criterion has already been met.
+"""
+function rgstep!(r::Renormalizer)
+    r.alg.stop(r.step, r.data) || error("stop criterion reached")
+    iterate(r, r.step)
+    return r.state, r.data
+end
+
+"""
+    run!(r::Renormalizer; verbosity=1)
+
+Run the RG flow to completion. Returns `(final_state, data)`.
+"""
+function run!(r::Renormalizer; verbosity = 1)
+    LoggingExtras.withlevel(; verbosity) do
+        @infov 1 "Starting simulation\n $(r.state)\n"
+        t = @elapsed for (_, data) in r
+            @infov 2 "Step $(r.step), data[end]: $(data[end])"
+        end
+        @infov 1 "Simulation finished\n $(stopping_info(r.alg.stop, r.step, r.data))\n Elapsed time: $(t)s\n Iterations: $(r.step)"
+    end
+    return r.state, r.data
+end
diff --git a/src/schemes/symmetric_looptnr.jl b/src/schemes/symmetric_looptnr.jl
index 2481cbae..691ce6f6 100644
--- a/src/schemes/symmetric_looptnr.jl
+++ b/src/schemes/symmetric_looptnr.jl
@@ -9,7 +9,7 @@ c4 & inversion symmetric Loop Optimization for Tensor Network Renormalization
 
 # Running the algorithm
     run!(::SLoopTNR, trscheme::TruncationStrategy,
-              criterion::TNRKit.stopcrit[, finalizer=default_Finalizer, finalize_beginning=true, oneloop=true,
+              criterion::TNRKit.stopcrit[, finalizer=default_Finalizer, oneloop=true,
               verbosity=1])
 
 `oneloop=true` will use disentangled tensors as a starting guess for the optimization.
@@ -194,17 +194,13 @@ end
 
 function run!(
         scheme::SLoopTNR, trscheme::TruncationStrategy,
-        criterion::TNRKit.stopcrit; finalizer = default_Finalizer, finalize_beginning = true, oneloop = true,
+        criterion::TNRKit.stopcrit; finalizer = default_Finalizer, oneloop = true,
         verbosity = 1
     )
     data = output_type(finalizer)[]
 
     LoggingExtras.withlevel(; verbosity) do
         @infov 1 "Starting simulation\n $(scheme)\n"
-
-        if finalize_beginning
-            push!(data, finalizer.f!(scheme))
-        end
         steps = 0
         crit = true
 
diff --git a/src/schemes/tnrscheme.jl b/src/schemes/tnrscheme.jl
index 00f5d365..8ae7177f 100644
--- a/src/schemes/tnrscheme.jl
+++ b/src/schemes/tnrscheme.jl
@@ -1,5 +1,4 @@
 # Extra code to make output type available
-function step! end
 function finalize! end
 
 """
@@ -10,7 +9,7 @@ Finalizer for TNR schemes
 # Constructors
     Finalizer(f!::Function, E::Type)
 
-A Finalizer holds a function `f!` that is to be applied to a TNR scheme after each step of the algorithm (and at the beginning if specified by `run!(;finalize_beginning=true)`, which is the default behavior).
+A Finalizer holds a function `f!` that is to be applied to a TNR scheme after each step of the algorithm.
 The type parameter `E` indicates the output type of `f!`, which is used to create an array of the correct type to hold the outputs.
 """
 struct Finalizer{E} # E is the output type of f
@@ -30,16 +29,14 @@ const ImpurityHOTRG_Finalizer = Finalizer(finalize!, Tuple{Float64, Float64, Flo
 # Finalization functions for the various TNR schemes
 abstract type TNRScheme{E, S} end
 
-function run!(scheme::TNRScheme, trscheme::TruncationStrategy, criterion::stopcrit, finalizer::Finalizer{E}; finalize_beginning = true, verbosity = 1) where {E}
+# TODO: This should be eventually replaced by `run!(renorm::Renormalizer; verbosity = 1)`.
+# Finalizers will be dropped at the end.
+function run!(scheme::TNRScheme, trscheme::TruncationStrategy, criterion::stopcrit, finalizer::Finalizer{E}; verbosity = 1) where {E}
     data = Vector{E}()
 
     LoggingExtras.withlevel(; verbosity) do
 
         @infov 1 "Starting simulation\n $(scheme)\n"
-        if finalize_beginning
-            push!(data, finalizer.f!(scheme))
-        end
-
         steps = 0
         crit = true
 
diff --git a/src/schemes/trg.jl b/src/schemes/trg.jl
index af1c9f43..de209042 100644
--- a/src/schemes/trg.jl
+++ b/src/schemes/trg.jl
@@ -1,23 +1,40 @@
 """
 $(TYPEDEF)
 
-Tensor Renormalization Group
+State holding a single tensor for a 1-site unit cell.
 
 # Constructors
     $(FUNCTIONNAME)(T)
 
-# Running the algorithm
-    run!(::TRG, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=default_Finalizer, finalize_beginning=true, verbosity=1])
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
+mutable struct OneSiteState{TT <: AbstractTensorMap{<:Any, <:Any, 2, 2}}
+    "central tensor"
+    T::TT
+
+    function OneSiteState(T::TT) where {TT <: AbstractTensorMap{<:Any, <:Any, 2, 2}}
+        return new{TT}(T)
+    end
+end
+
+"""
+$(TYPEDEF)
+
+Tensor Renormalization Group algorithm.
 
 Each step rescales the lattice by a (linear) factor of √2,
-and rotate the lattice by 45 degrees in counter clockwise direction.
+and rotates the lattice by 45 degrees in counter clockwise direction.
 The elementary modular parameter `τ₀ ↦ (τ₀ - 1) / (τ₀ + 1)`.
 
-!!! info "verbosity levels"
-    - 0: No output
-    - 1: Print information at start and end of the algorithm
-    - 2: Print information at each step
-    
+# Constructors
+
+All parameters are passed as keyword arguments with sensible defaults:
+
+    $(FUNCTIONNAME)(; trunc=truncrank(16), stop=maxiter(20))
+
 # Fields
 
 $(TYPEDFIELDS)
@@ -25,25 +42,33 @@ $(TYPEDFIELDS)
 # References
 * [Levin & Nave Phys. Rev. Letters 99(12) (2007)](@cite levin2007)
 """
-mutable struct TRG{E, S, TT <: AbstractTensorMap{E, S, 2, 2}} <: TNRScheme{E, S}
-    "central tensor"
-    T::TT
+Base.@kwdef struct TRG <: TNRAlgorithm
+    "Truncation strategy for SVD steps"
+    trunc::TruncationStrategy = truncrank(16)
+    "Stopping criterion"
+    stop::stopcrit = maxiter(20)
+end
 
-    function TRG(T::TT) where {E, S, TT <: AbstractTensorMap{E, S, 2, 2}}
-        return new{E, S, TT}(T)
-    end
+function Renormalizer(alg::TRG, T::TT) where {TT}
+    state = OneSiteState(T)
+    return Renormalizer(alg, state, scalartype(T)[])
+end
+
+function Renormalizer(alg::TRG, state::OneSiteState)
+    return Renormalizer(alg, state, scalartype(state.T)[])
 end
 
-function step!(scheme::TRG, trunc::TruncationStrategy)
-    A, B = SVD12(scheme.T, trunc)
-    Tp = transpose(scheme.T, ((2, 4), (1, 3)))
-    C, D = SVD12(Tp, trunc)
-    @plansor scheme.T[-1 -2; -3 -4] := D[-2; 1 2] * B[-1; 4 1] * C[4 3; -3] * A[3 2; -4]
-    return scheme
+function step!(state::OneSiteState, alg::TRG)
+    A, B = SVD12(state.T, alg.trunc)
+    Tp = transpose(state.T, ((2, 4), (1, 3)))
+    C, D = SVD12(Tp, alg.trunc)
+    @plansor state.T[-1 -2; -3 -4] := D[-2; 1 2] * B[-1; 4 1] * C[4 3; -3] * A[3 2; -4]
+    return state
 end
 
-function Base.show(io::IO, scheme::TRG)
-    println(io, "TRG - Tensor Renormalization Group")
-    println(io, "  * T: $(summary(scheme.T))")
+function Base.show(io::IO, alg::TRG)
+    println(io, "TRG")
+    println(io, "  * truncation: $(alg.trunc)")
+    println(io, "  * stop: $(alg.stop)")
     return nothing
 end
diff --git a/src/schemes/ttnr.jl b/src/schemes/ttnr.jl
index 2d2d2a62..f2110e05 100644
--- a/src/schemes/ttnr.jl
+++ b/src/schemes/ttnr.jl
@@ -63,7 +63,7 @@ Thermal Tensor Network Renormalization on a square-lattice
 
 # Running the algorithm
     run!(::ThermalTNR, A::AbstractMatrix{<:AbstractTensorMap}, trunc::TruncationStrategy, criterion::stopcrit[
-              , finalizer=default_Finalizer, finalize_beginning=true, verbosity=1])
+              , finalizer=default_Finalizer, verbosity=1])
 
 # Fields
 
@@ -255,17 +255,12 @@ end
 
 function run!(
         scheme::ThermalTNR, layer::TNO, trscheme::TruncationStrategy,
-        criterion::stopcrit, finalizer::Finalizer{E};
-        finalize_beginning = true, verbosity = 1
+        criterion::stopcrit, finalizer::Finalizer{E}; verbosity = 1
     ) where {E}
     data = Vector{E}()
 
     LoggingExtras.withlevel(; verbosity) do
         @infov 1 "Starting simulation\n $(scheme)\n"
-        if finalize_beginning
-            push!(data, finalizer.f!(scheme))
-        end
-
         steps = 0
         crit = true
 
diff --git a/src/utility/finalize.jl b/src/utility/finalize.jl
index ec271ed1..54dfdb01 100644
--- a/src/utility/finalize.jl
+++ b/src/utility/finalize.jl
@@ -1,4 +1,17 @@
-const simple_scheme = Union{TRG, ATRG, HOTRG}
+# ========================================================
+# normalization after an RG step for different TNR schemes
+# ========================================================
+
+# 1x1 unitcell finalize
+function finalize!(state::OneSiteState)
+    n = norm(@tensor state.T[1 2; 2 1])
+    state.T /= n
+    return n
+end
+
+# Below: for use with old `TNRScheme` interface.
+
+const simple_scheme = Union{ATRG, HOTRG}
 
 # 1x1 unitcell finalize
 function finalize!(scheme::simple_scheme)
@@ -134,6 +147,10 @@ function finalize_phase23!(scheme::CorrelationHOTRG)
     return n, n_imp, n_imp
 end
 
+# =====================================================================
+# TODO: remove the following once Renormalizer interface is finished
+# =====================================================================
+
 # cft data finalize
 function finalize_cftdata!(scheme::TNRScheme)
     finalize!(scheme)
diff --git a/src/utility/free_energy.jl b/src/utility/free_energy.jl
index 65604d04..43810b89 100644
--- a/src/utility/free_energy.jl
+++ b/src/utility/free_energy.jl
@@ -12,7 +12,7 @@ and computes the free energy.
 """
 function free_energy(data, β; scalefactor = 2.0, initial_size = 1.0)
     lnz = 0.0
-    x = 1.0 - log(initial_size) / log(scalefactor)
+    x = -log(initial_size) / log(scalefactor)
     for (i, z) in enumerate(data)
         lnz += log(z) * scalefactor^(x - i)
     end
diff --git a/src/utility/gs_degeneracy.jl b/src/utility/gs_degeneracy.jl
index cffd000f..2be2c2b0 100644
--- a/src/utility/gs_degeneracy.jl
+++ b/src/utility/gs_degeneracy.jl
@@ -4,147 +4,77 @@
 Calculates the Ground State Degeneracy (GSD) from the fixed-point tensor of a TNRScheme,
 using the eigenvalues of the transfer matrix. The GSD is the exponential of the Shannon entropy.
 """
-function ground_state_degeneracy(scheme::TNRScheme{E}, unitcell::Int = 1) where {E}
-    # Construct contraction indices
-    indices = Vector{NTuple{4, Int}}(undef, unitcell)
-    for i in 1:unitcell
-        indices[i] = (i, -i, -(i + unitcell), i + 1)
-    end
-    indices[end] = (unitcell, -unitcell, -(unitcell + unitcell), 1)
-
-    # Contract tensors
-    Ts = fill(scheme.T, unitcell)
-    T = ncon(Ts, indices)
-
-    # Construct static tuple indices
-    outinds = ntuple(i -> i, unitcell)
-    ininds = ntuple(i -> unitcell + i, unitcell)
-
-    T = permute(T, (outinds, ininds))
-
-    # Compute normalized eigenvalues
-    D, _ = eig_full(T)
-    D = D / tr(D)
-    vals = filter(!iszero, abs.(D.data))
-    # Shannon entropy (stable + efficient)
-    S = 0.0
-    for v in vals
-        ev = abs(v)
-        if ev > 0
-            S -= ev * log(ev)
-        end
-    end
-
-    return exp(S)
+function ground_state_degeneracy(T::AbstractTensorMap, unitcell::Int = 1)
+    tm = _row_transfer_matrix(T, unitcell)
+    return _ground_state_degeneracy(tm)
 end
-function ground_state_degeneracy(scheme::BTRG{E}; unitcell::Int = 1) where {E}
-    indices = Vector{NTuple{4, Int}}(undef, unitcell)
-    for i in 1:unitcell
-        indices[i] = (i, -i, -(i + unitcell), i + 1)
-    end
-    indices[end] = (unitcell, -unitcell, -(unitcell + unitcell), 1)
-
+function ground_state_degeneracy(TA::AbstractTensorMap, TB::AbstractTensorMap)
+    # 2-column transfer matrix
+    @tensor tm[-1 -2; -3 -4] := TA[-1 1; 3 2] * TB[2 6; 4 -3] *
+        TB[-2 3; 1 5] * TA[5 4; 6 -4]
+    return _ground_state_degeneracy(tm)
+end
+ground_state_degeneracy(scheme::TNRScheme; unitcell::Int = 1) = ground_state_degeneracy(scheme.T, unitcell)
+function ground_state_degeneracy(scheme::BTRG; unitcell::Int = 1)
     @tensor T_unit[-1 -2; -3 -4] := scheme.T[1 2; -3 -4] * scheme.S1[-2; 2] *
         scheme.S2[-1; 1]
-    T = ncon(fill(T_unit, unitcell), indices)
-
-    # Construct static tuple indices
-    outinds = ntuple(i -> i, unitcell)
-    ininds = ntuple(i -> unitcell + i, unitcell)
-
-    T = permute(T, (outinds, ininds))
-    D, _ = eig_full(T)
-    D = D / tr(D)
-    vals = filter(!iszero, abs.(D.data))
-    # Shannon entropy (stable + efficient)
-    S = 0.0
-    for v in vals
-        ev = abs(v)
-        if ev > 0
-            S -= ev * log(ev)
-        end
-    end
-
-    return exp(S)
+    return ground_state_degeneracy(T_unit, unitcell)
 end
-function ground_state_degeneracy(scheme::LoopTNR{E}) where {E}
-    norm_const = area_term(scheme.TA, scheme.TB)
-    T1 = scheme.TA / abs(norm_const)^(1 / 4)
-    T2 = scheme.TB / abs(norm_const)^(1 / 4)
 
-    @tensor T_unit[-1 -2; -3 -4] := T1[-1 1; 3 2] * T2[2 6; 4 -3] *
-        T2[-2 3; 1 5] * T1[5 4; 6 -4]
-
-    D, _ = eig_full(T_unit)
+ground_state_degeneracy(scheme::LoopTNR) = ground_state_degeneracy(scheme.TA, scheme.TB)
+# helper function
+function _ground_state_degeneracy(tm::AbstractTensorMap{E, S, N, N}) where {E, S, N}
+    D, _ = eig_full(tm)
     D = D / tr(D)
-    vals = filter(!iszero, abs.(D.data))
-    # Shannon entropy (stable + efficient)
-    S = 0.0
-    for v in vals
-        ev = abs(v)
-        if ev > 0
-            S -= ev * log(ev)
-        end
-    end
-
-    return exp(S)
+    evs = filter(!iszero, abs.(D.data))
+    entropy = -sum(evs .* log.(evs))
+    return exp(entropy)
 end
 
 """
 $(SIGNATURES)
     
-Calculates the Gu-Wen ratio X1 and X2 from the fixed-point tensor of a TNRScheme.
+Calculates the Gu-Wen ratio X1 and X2 from the fixed-point tensor(s).
 The Gu-Wen ratios are related to the Ground state Degeneracy and the the scaling dimensions. See references.
 
 # References
 * [Zheng-Cheng Gu & Xiao-Gang Wen. PhysRevB.80.155131](@cite gu2009)
 * [Satoshi Morita et al. arxiv:2512.03395](@cite morita2025)
 """
-function gu_wen_ratio(scheme::TNRScheme{E}) where {E}
-    T_unit = scheme.T
-
-    one_norm = norm(@tensor T_unit[1 2; 2 1])
-    two_norm_X1 = norm(@tensor T_unit[1 2; 2 3] * T_unit[3 4; 4 1])
-    two_norm_X2 = norm(@tensor T_unit[1 2; 3 4] * T_unit[4 3; 2 1])
-
+function gu_wen_ratio(T::AbstractTensorMap{E, S, 2, 2}) where {E, S}
+    one_norm = norm(@tensor T[1 2; 2 1])
+    two_norm_X1 = norm(@tensor T[1 2; 2 3] * T[3 4; 4 1])
+    two_norm_X2 = norm(@tensor T[1 2; 3 4] * T[4 3; 2 1])
     X1 = (one_norm^2) / (two_norm_X1)
     X2 = (one_norm^2) / (two_norm_X2)
     return X1, X2
 end
-function gu_wen_ratio(scheme::BTRG{E}) where {E}
-    @tensor T_unit[-1 -2; -3 -4] := scheme.T[1 2; -3 -4] * scheme.S1[-2; 2] *
-        scheme.S2[-1; 1]
-
-    one_norm = norm(@tensor T_unit[1 2; 2 1])
-    two_norm_X1 = norm(@tensor T_unit[1 2; 2 3] * T_unit[3 4; 4 1])
-    two_norm_X2 = norm(@tensor T_unit[1 2; 3 4] * T_unit[4 3; 2 1])
-
-    X1 = (one_norm^2) / (two_norm_X1)
-    X2 = (one_norm^2) / (two_norm_X2)
-    return X1, X2
-end
-function gu_wen_ratio(scheme::LoopTNR{E}) where {E}
-    T1 = scheme.TA
-    T2 = scheme.TB
+function gu_wen_ratio(
+        TA::AbstractTensorMap{E, S, 2, 2}, TB::AbstractTensorMap{E, S, 2, 2}
+    ) where {E, S}
     one_norm = norm(
-        @tensor opt = true T1[1 2; 3 4] * T2[4 5; 6 1] *
-            T2[7 3; 2 8] * T1[8 6; 5 7]
+        @tensor opt = true TA[1 2; 3 4] * TB[4 5; 6 1] *
+            TB[7 3; 2 8] * TA[8 6; 5 7]
     )
-
     two_norm_X1 = norm(
-        @tensor opt = true T1[1 2; 3 4] * T2[4 5; 6 7] *
-            T1[7 8; 9 10] * T2[10 11; 12 1] *
-            T2[13 3; 2 14] * T1[14 6; 5 15] * T2[15 9; 8 16] * T1[16 12; 11 13]
+        @tensor opt = true TA[1 2; 3 4] * TB[4 5; 6 7] *
+            TA[7 8; 9 10] * TB[10 11; 12 1] *
+            TB[13 3; 2 14] * TA[14 6; 5 15] * TB[15 9; 8 16] * TA[16 12; 11 13]
     )
-
     two_norm_X2 = norm(
-        @tensor opt = true T1[1 2; 3 4] * T2[4 5; 6 7] *
-            T1[7 8; 9 10] * T2[10 11; 12 1] *
-            T2[13 9; 2 14] * T1[14 12; 5 15] *
-            T2[15 3; 8 16] * T1[16 6; 11 13]
+        @tensor opt = true TA[1 2; 3 4] * TB[4 5; 6 7] *
+            TA[7 8; 9 10] * TB[10 11; 12 1] *
+            TB[13 9; 2 14] * TA[14 12; 5 15] *
+            TB[15 3; 8 16] * TA[16 6; 11 13]
     )
-
     X1 = (one_norm^2) / (two_norm_X1)
     X2 = (one_norm^2) / (two_norm_X2)
     return X1, X2
 end
+gu_wen_ratio(scheme::TNRScheme) = gu_wen_ratio(scheme.T)
+function gu_wen_ratio(scheme::BTRG)
+    @tensor T_unit[-1 -2; -3 -4] := scheme.T[1 2; -3 -4] * scheme.S1[-2; 2] *
+        scheme.S2[-1; 1]
+    return gu_wen_ratio(T_unit)
+end
+gu_wen_ratio(scheme::LoopTNR) = gu_wen_ratio(scheme.TA, scheme.TB)
diff --git a/test/fermions/fermions.jl b/test/fermions/fermions.jl
index f434c503..cbcb3762 100644
--- a/test/fermions/fermions.jl
+++ b/test/fermions/fermions.jl
@@ -7,8 +7,9 @@ T = gross_neveu_start(0, 0, 0)
 
 # === TRG ===
 @testset "TRG - Gross-Neveu Model" begin
-    scheme = TRG(T)
-    data = run!(scheme, truncrank(16), maxiter(25))
+    alg = TRG(; trunc = truncrank(16), stop = maxiter(25))
+    renorm = Renormalizer(alg, T)
+    _, data = run!(renorm; verbosity = 0)
     @test free_energy(data, 1.0) ≈ f_bench rtol = 1.0e-3
 end
 
@@ -56,7 +57,14 @@ end
     data_c4vCTM = run!(c4vCTM(T_flipped_C4v), truncrank(8), maxiter(10))
     free_energy_c4vCTM = -data_c4vCTM / β
 
-    schemes = [TRG, BTRG, HOTRG, ATRG, LoopTNR]
+    # TRG with new iterable interface
+    alg = TRG(; trunc = truncrank(8), stop = maxiter(10))
+    renorm = Renormalizer(alg, T_flipped_C4v)
+    _, data = run!(renorm; verbosity = 0)
+    @test free_energy_c4vCTM ≈ free_energy(data, β; scalefactor = 2.0) rtol = 1.0e-9
+
+    # old interface
+    schemes = [BTRG, HOTRG, ATRG, LoopTNR]
     for scheme in schemes
         data = run!(scheme(T_flipped_C4v), truncrank(8), maxiter(10))
         scalefactor = scheme ∈ [HOTRG, ATRG] ? 4.0 : 2.0
diff --git a/test/models/models.jl b/test/models/models.jl
index d4ccbbd5..b0ed3e29 100644
--- a/test/models/models.jl
+++ b/test/models/models.jl
@@ -44,8 +44,9 @@ model_temp_answer_string_3d = [
 
 for (model, temp, answer, description) in model_temp_answer_string_2d
     @testset "$(description)" begin
-        scheme = TRG(model)
-        data = run!(scheme, truncrank(16), maxiter(25))
+        alg = TRG(; trunc = truncrank(16), stop = maxiter(25))
+        renorm = Renormalizer(alg, model)
+        _, data = run!(renorm; verbosity = 0)
         @test free_energy(data, temp) ≈ answer rtol = 1.0e-3
     end
 end
@@ -201,7 +202,7 @@ end
     scheme = LoopTNR(T)
     elt = scalartype(T)
     finalizer = Finalizer(cc_finalize!, Tuple{elt, complex(elt), elt})
-    data = run!(scheme, truncrank(16), maxiter(16), finalizer; finalize_beginning = false)
+    data = run!(scheme, truncrank(16), maxiter(16), finalizer)
     @test last(data)[3] ≈ 0.5 atol = 1.0e-2
 end
 
@@ -213,9 +214,9 @@ end
     T = vertical_stack_exp(T, nfold, trunc_stack)
 
     @info "CFT data from TRG"
-    scheme = TRG(T)
-    run!(scheme, truncrank(24), maxiter(8))
-    cft = CFTData(scheme; shape = [1, 1, 0])
+    alg = TRG(; trunc = truncrank(24), stop = maxiter(8))
+    state, _ = run!(Renormalizer(alg, T); verbosity = 0)
+    cft = CFTData(state.T; shape = [1, 1, 0])
     c = cft.central_charge
     sd = cft.scaling_dimensions
     d_1 = real(sd[(:NS, FermionParity(0))][2])
@@ -233,7 +234,7 @@ end
     scheme = LoopTNR(T)
     elt = scalartype(T)
     finalizer = Finalizer(cc_finalize!, Tuple{elt, complex(elt), elt})
-    data = run!(scheme, truncrank(16), maxiter(8), finalizer; finalize_beginning = false)
+    data = run!(scheme, truncrank(16), maxiter(8), finalizer)
     @test last(data)[3] ≈ 0.5 atol = 1.0e-2
     for shape in ([√2, 2√2, 0], [1, 4, 1], [1, 8, 1], [4 / √10, 2√10, 2 / √10])
         cft = CFTData(scheme; shape = shape)
diff --git a/test/schemes/schemes.jl b/test/schemes/schemes.jl
index b428a0a6..cf082ad2 100644
--- a/test/schemes/schemes.jl
+++ b/test/schemes/schemes.jl
@@ -17,19 +17,9 @@ const T_aniso = classical_ising(Z2Irrep, βc_aniso; Jx = Jx_aniso, Jy = Jy_aniso
 const f_aniso_exact = f_onsager_anisotropic(βc_aniso, Jx_aniso, Jy_aniso)
 const τ_aniso_exact = sinh(2 * βc_aniso * Jx_aniso)
 
-function cft_finalize!(scheme)
-    finalize!(scheme)
-    return CFTData(scheme)
-end
-
 """
 Normalize the tensor, return the normalization factor and elementary modular parameter
 """
-function tau_finalize!(scheme::TRG)
-    n = finalize!(scheme)
-    τ0, c = extract_tau_and_c(scheme.T; fast = false)
-    return (n, τ0)
-end
 function tau_finalize!(scheme::LoopTNR)
     n = finalize!(scheme)
     τ0, c = extract_tau_and_c(scheme.TA, scheme.TB; fast = false)
@@ -40,27 +30,28 @@ end
 @testset "TRG - Anisotropic Ising Model" begin
     @info "Anisotropy: Jx = $(Jx_aniso), Jy = $(Jy_aniso)"
     @info "TRG anisotropic ising free energy"
-    scheme = TRG(T_aniso)
-    elt = scalartype(T_aniso)
-    finalizer = Finalizer(tau_finalize!, Tuple{elt, complex(elt)})
-    data = run!(scheme, truncrank(24), maxiter(25), finalizer)
-
-    ns = map(Base.Fix2(getindex, 1), data)
-    @test free_energy(ns, βc_aniso) ≈ f_aniso_exact rtol = 2.0e-6
+    alg = TRG(; trunc = truncrank(24), stop = maxiter(25))
+    renorm = Renormalizer(alg, T_aniso)
+    T_step10 = nothing
+    τs = complex(scalartype(T_aniso))[]
+    for (state, _) in renorm
+        τ = first(extract_tau_and_c(state.T; fast = false))
+        @info "* Step $(renorm.step): τ = $τ."
+        push!(τs, τ)
+        (renorm.step == 10) && (T_step10 = state.T)
+    end
+    @test free_energy(renorm.data, βc_aniso) ≈ f_aniso_exact rtol = 2.0e-6
 
     @info "TRG τ → (τ - 1) / (τ + 1)"
     f_trg(τ) = (τ - 1) / (τ + 1)
-    τs = map(Base.Fix2(getindex, 2), data)
     for n in 5:7
         @test τs[n + 1] ≈ f_trg(τs[n]) rtol = 5.0e-2
         @info "* verified for step $(n - 1) → $n"
     end
 
     @info "TRG anisotropic ising CFT data — shape [1, 1, 0]"
-    scheme = TRG(T_aniso)
-    run!(scheme, truncrank(24), maxiter(10))
     # use fast tau algorithm below
-    cft = CFTData(scheme; shape = [1, 1, 0])
+    cft = CFTData(T_step10; shape = [1, 1, 0])
     sd_all = real(cft.scaling_dimensions)
     cft_sorted = sort(sd_all[2:end]; by = abs)
 
@@ -69,23 +60,19 @@ end
     @info "Obtained scaling dimensions: Δ₁ = $(cft_sorted[1]), Δ₂ = $(cft_sorted[2])"
 
     @info "TRG anisotropic ising ground state degeneracy"
-    T1 = classical_ising(βc_aniso - 0.01; Jx = Jx_aniso, Jy = Jy_aniso)
-    scheme = TRG(T1)
-    run!(scheme, truncrank(16), maxiter(20))
-    gsd = ground_state_degeneracy(scheme)
-    X1, X2 = gu_wen_ratio(scheme)
-    @test gsd ≈ 1 rtol = 1.0e-2
-    @test X1 ≈ 1.0 rtol = 1.0e-2
-    @test X2 ≈ 1.0 rtol = 1.0e-2
-
-    T2 = classical_ising(βc_aniso + 0.01; Jx = Jx_aniso, Jy = Jy_aniso)
-    scheme = TRG(T2)
-    run!(scheme, truncrank(16), maxiter(20))
-    gsd = ground_state_degeneracy(scheme)
-    X1, X2 = gu_wen_ratio(scheme)
-    @test gsd ≈ 2 rtol = 1.0e-2
-    @test X1 ≈ 2.0 rtol = 1.0e-2
-    @test X2 ≈ 2.0 rtol = 1.0e-2
+    alg = TRG(; trunc = truncrank(16), stop = maxiter(20))
+    for (Δβ, gsd0, X10, X20) in [
+            (-0.01, 1, 1.0, 1.0), (0.01, 2, 2.0, 2.0),
+        ]
+        T1 = classical_ising(βc_aniso + Δβ; Jx = Jx_aniso, Jy = Jy_aniso)
+        renorm = Renormalizer(alg, T1)
+        state, = run!(renorm; verbosity = 0)
+        gsd = ground_state_degeneracy(state.T)
+        X1, X2 = gu_wen_ratio(state.T)
+        @test gsd ≈ gsd0 rtol = 1.0e-2
+        @test X1 ≈ X10 rtol = 1.0e-2
+        @test X2 ≈ X20 rtol = 1.0e-2
+    end
 end
 
 # BTRG
@@ -106,23 +93,18 @@ end
     @test cft[2] ≈ ising_cft_exact[2] rtol = 2.0e-2
 
     @info "BTRG ising ground state degeneracy"
-    T1 = classical_ising(ising_βc - 0.01)
-    scheme = BTRG(T1)
-    run!(scheme, truncrank(16), maxiter(20))
-    gsd = ground_state_degeneracy(scheme)
-    X1, X2 = gu_wen_ratio(scheme)
-    @test gsd ≈ 1 rtol = 1.0e-2
-    @test X1 ≈ 1.0 rtol = 1.0e-2
-    @test X2 ≈ 1.0 rtol = 1.0e-2
-
-    T2 = classical_ising(ising_βc + 0.01)
-    scheme = BTRG(T2)
-    run!(scheme, truncrank(16), maxiter(20))
-    gsd = ground_state_degeneracy(scheme)
-    X1, X2 = gu_wen_ratio(scheme)
-    @test gsd ≈ 2 rtol = 1.0e-2
-    @test X1 ≈ 2.0 rtol = 1.0e-2
-    @test X2 ≈ 2.0 rtol = 1.0e-2
+    for (Δβ, gsd0, X10, X20) in [
+            (-0.01, 1, 1.0, 1.0), (0.01, 2, 2.0, 2.0),
+        ]
+        T1 = classical_ising(ising_βc + Δβ)
+        scheme = BTRG(T1)
+        run!(scheme, truncrank(16), maxiter(20))
+        gsd = ground_state_degeneracy(scheme)
+        X1, X2 = gu_wen_ratio(scheme)
+        @test gsd ≈ gsd0 rtol = 1.0e-2
+        @test X1 ≈ X10 rtol = 1.0e-2
+        @test X2 ≈ X20 rtol = 1.0e-2
+    end
 end
 
 # HOTRG
@@ -143,23 +125,18 @@ end
     @test cft[2] ≈ ising_cft_exact[2] rtol = 1.0e-2
 
     @info "HOTRG ising ground state degeneracy"
-    T1 = classical_ising(ising_βc - 0.01)
-    scheme = HOTRG(T1)
-    run!(scheme, truncrank(12), maxiter(20))
-    gsd = ground_state_degeneracy(scheme)
-    X1, X2 = gu_wen_ratio(scheme)
-    @test gsd ≈ 1 rtol = 1.0e-2
-    @test X1 ≈ 1.0 rtol = 1.0e-2
-    @test X2 ≈ 1.0 rtol = 1.0e-2
-
-    T2 = classical_ising(ising_βc + 0.01)
-    scheme = HOTRG(T2)
-    run!(scheme, truncrank(12), maxiter(20))
-    gsd = ground_state_degeneracy(scheme)
-    X1, X2 = gu_wen_ratio(scheme)
-    @test gsd ≈ 2 rtol = 1.0e-2
-    @test X1 ≈ 2.0 rtol = 1.0e-2
-    @test X2 ≈ 2.0 rtol = 1.0e-2
+    for (Δβ, gsd0, X10, X20) in [
+            (-0.01, 1, 1.0, 1.0), (0.01, 2, 2.0, 2.0),
+        ]
+        T1 = classical_ising(ising_βc + Δβ)
+        scheme = HOTRG(T1)
+        run!(scheme, truncrank(12), maxiter(20))
+        gsd = ground_state_degeneracy(scheme)
+        X1, X2 = gu_wen_ratio(scheme)
+        @test gsd ≈ gsd0 rtol = 1.0e-2
+        @test X1 ≈ X10 rtol = 1.0e-2
+        @test X2 ≈ X20 rtol = 1.0e-2
+    end
 end
 
 # ATRG
@@ -180,23 +157,18 @@ end
     @test cft[2] ≈ ising_cft_exact[2] rtol = 1.0e-2
 
     @info "ATRG ising ground state degeneracy"
-    T1 = classical_ising(ising_βc - 0.01)
-    scheme = ATRG(T1)
-    run!(scheme, truncrank(16), maxiter(20))
-    gsd = ground_state_degeneracy(scheme)
-    X1, X2 = gu_wen_ratio(scheme)
-    @test gsd ≈ 1 rtol = 1.0e-2
-    @test X1 ≈ 1.0 rtol = 1.0e-2
-    @test X2 ≈ 1.0 rtol = 1.0e-2
-
-    T2 = classical_ising(ising_βc + 0.01)
-    scheme = ATRG(T2)
-    run!(scheme, truncrank(16), maxiter(20))
-    gsd = ground_state_degeneracy(scheme)
-    X1, X2 = gu_wen_ratio(scheme)
-    @test gsd ≈ 2 rtol = 1.0e-2
-    @test X1 ≈ 2.0 rtol = 1.0e-2
-    @test X2 ≈ 2.0 rtol = 1.0e-2
+    for (Δβ, gsd0, X10, X20) in [
+            (-0.01, 1, 1.0, 1.0), (0.01, 2, 2.0, 2.0),
+        ]
+        T1 = classical_ising(ising_βc + Δβ)
+        scheme = ATRG(T1)
+        run!(scheme, truncrank(16), maxiter(20))
+        gsd = ground_state_degeneracy(scheme)
+        X1, X2 = gu_wen_ratio(scheme)
+        @test gsd ≈ gsd0 rtol = 1.0e-2
+        @test X1 ≈ X10 rtol = 1.0e-2
+        @test X2 ≈ X20 rtol = 1.0e-2
+    end
 end
 
 # LoopTNR
@@ -222,13 +194,12 @@ end
     τs = map(Base.Fix2(getindex, 2), data)
     for n in 5:8
         @test τs[n + 1] ≈ f_looptnr(τs[n]) rtol = 2.0e-2
-        @info "* verified for step $(n - 1) → $n"
+        @info "* verified for step $n → $(n + 1)"
     end
 
     @info "Theory value of τ for anisotropic Ising"
     for n in (4, 8, 12)
-        # n + 1 due to finalizing the initial tensor
-        τ = τs[n + 1]
+        τ = τs[n]
         @test real(τ) ≈ 0 atol = 1.0e-3
         @test imag(τ) ≈ τ_aniso_exact rtol = 2.0e-3
         @info "* τ = $τ ≈ $(τ_aniso_exact)im verified for step $n"
@@ -282,23 +253,18 @@ end
     end
 
     @info "LoopTNR anisotropic ising ground state degeneracy"
-    T1 = classical_ising(βc_aniso - 0.01; Jx = Jx_aniso, Jy = Jy_aniso)
-    scheme = LoopTNR(T1)
-    run!(scheme, truncrank(12), maxiter(20))
-    gsd = ground_state_degeneracy(scheme)
-    X1, X2 = gu_wen_ratio(scheme)
-    @test gsd ≈ 1 rtol = 1.0e-2
-    @test X1 ≈ 1.0 rtol = 1.0e-2
-    @test X2 ≈ 1.0 rtol = 1.0e-2
-
-    T2 = classical_ising(βc_aniso + 0.01; Jx = Jx_aniso, Jy = Jy_aniso)
-    scheme = LoopTNR(T2)
-    run!(scheme, truncrank(12), maxiter(20))
-    gsd = ground_state_degeneracy(scheme)
-    X1, X2 = gu_wen_ratio(scheme)
-    @test gsd ≈ 2 rtol = 1.0e-2
-    @test X1 ≈ 2.0 rtol = 1.0e-2
-    @test X2 ≈ 2.0 rtol = 1.0e-2
+    for (Δβ, gsd0, X10, X20) in [
+            (-0.01, 1, 1.0, 1.0), (0.01, 2, 2.0, 2.0),
+        ]
+        T1 = classical_ising(βc_aniso + Δβ; Jx = Jx_aniso, Jy = Jy_aniso)
+        scheme = LoopTNR(T1)
+        run!(scheme, truncrank(12), maxiter(20))
+        gsd = ground_state_degeneracy(scheme)
+        X1, X2 = gu_wen_ratio(scheme)
+        @test gsd ≈ gsd0 rtol = 1.0e-2
+        @test X1 ≈ X10 rtol = 1.0e-2
+        @test X2 ≈ X20 rtol = 1.0e-2
+    end
 end
 
 @testset "LoopTNR - Ising Model - Dense Solver - NNR" begin
@@ -570,7 +536,7 @@ end
 
     data = run!(scheme, truncrank(16), maxiter(25))
 
-    @test free_energy(getindex.(data, 1), ising_βc; scalefactor = 4.0, initial_size = 4.0) ≈ f_onsager rtol = 1.0e-4
+    @test free_energy(getindex.(data, 1), ising_βc; scalefactor = 4.0) ≈ f_onsager rtol = 1.0e-4
 end
 
 @testset "Correlation HOTRG - Magnetisation Correlation" begin
@@ -728,8 +694,7 @@ function _thermal_zn2_gu_wen_x1(β, Lz; χttnr = 12, χbtrg = 16, btrg_steps = 1
     @tensor effective_tensor[-1 -2 -3 -4] := scheme.T[1, 1][p p; -1 -2 -3 -4]
     btrg = BTRG(permute(effective_tensor, ((1, 2), (3, 4))))
     ratios = run!(
-        btrg, truncrank(χbtrg), maxiter(btrg_steps), guwenratio_Finalizer;
-        finalize_beginning = false, verbosity = 0,
+        btrg, truncrank(χbtrg), maxiter(btrg_steps), guwenratio_Finalizer; verbosity = 0,
     )
 
     x1, _ = last(ratios)