From 46a86d1addff52fea18231604da30768c328c99a Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Mon, 20 Jul 2026 14:35:59 +0200
Subject: [PATCH 01/19] feat: add post expansion algorithms

---
 src/MPSKit.jl                             |  5 ++
 src/algorithms/post_expand/dmrg3s.jl      | 72 +++++++++++++++++++++++
 src/algorithms/post_expand/post_expand.jl | 11 ++++
 3 files changed, 88 insertions(+)
 create mode 100644 src/algorithms/post_expand/dmrg3s.jl
 create mode 100644 src/algorithms/post_expand/post_expand.jl

diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index 2ef8cf900..3240f9a25 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -38,6 +38,8 @@ export time_evolve, timestep, timestep!, make_time_mpo
 export TDVP, TDVP2, WI, WII, TaylorCluster
 export changebonds, changebonds!
 export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand
+export post_expand!, NoExpand
+export NoiseSchedule, ExponentialDecay, Warmup, DMRG3S
 export propagator
 export DynamicalDMRG, NaiveInvert, Jeckelmann
 export exact_diagonalization, fidelity_susceptibility
@@ -160,6 +162,9 @@ include("algorithms/changebonds/svdcut.jl")
 include("algorithms/changebonds/randexpand.jl")
 include("algorithms/changebonds/sketchedexpand.jl")
 
+include("algorithms/post_expand/post_expand.jl")
+include("algorithms/post_expand/dmrg3s.jl")
+
 include("algorithms/timestep/tdvp.jl")
 include("algorithms/timestep/taylorcluster.jl")
 include("algorithms/timestep/wii.jl")
diff --git a/src/algorithms/post_expand/dmrg3s.jl b/src/algorithms/post_expand/dmrg3s.jl
new file mode 100644
index 000000000..9e613c814
--- /dev/null
+++ b/src/algorithms/post_expand/dmrg3s.jl
@@ -0,0 +1,72 @@
+abstract type NoiseSchedule end
+
+struct ExponentialDecay <: NoiseSchedule
+    decay_rate::Float64
+end
+
+struct Warmup <: NoiseSchedule
+    iters::Int
+end
+
+(s::ExponentialDecay)(noise, iter, ϵ) = noise * s.decay_rate^iter
+(s::Warmup)(noise, iter, ϵ) = iter ≤ s.iters ? noise : zero(noise)
+
+struct DMRG3S{N, S <: NoiseSchedule} <: Algorithm
+    noise::N
+    schedule::S
+end
+
+function _update_post_expand(alg::DMRG3S, iter, ϵ)
+    noise = alg.schedule(alg.noise, iter, ϵ)
+    return iszero(noise) ? NoExpand() : DMRG3S(noise, alg.schedule)
+end
+
+function _get_combiner(::Type{T}, V1, V2) where {T}
+    Vf = fuse(V1 ⊗ V2)
+    return isomorphism(T, (V1 ⊗ V2) ← Vf), Vf
+end
+
+function post_expand!(pos::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::DMRG3S, envs, AC, alg_gauge; normalize=true)
+    El = leftenv(envs, pos, ψ)
+    Hi = H[pos]
+    α = alg.noise
+    T = promote_type(scalartype(ψ), scalartype(Hi))
+    V = right_virtualspace(AC)
+    combiner, Vpert = _get_combiner(T, V, right_virtualspace(Hi))
+    
+    @plansor pert[-1 -2; -3] := α * El[-1 1; 2] * AC[2 3; 4] * Hi[1 -2; 3 5] * combiner[4, 5; -3]    
+
+    AC_expanded = catdomain(AC, pert)
+
+    AL, C = left_gauge(AC_expanded, alg_gauge)
+    B = _transpose_tail(ψ.AR[pos+1])
+    AR = _transpose_front(catcodomain(B, zeros(T, Vpert ← domain(B))))
+
+    normalize && normalize!(C)
+    ψ.AC[pos] = (AL, C)
+    ψ.AC[pos + 1] = (C, AR)
+    return ψ
+end
+
+function post_expand!(pos::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::DMRG3S, envs, AC, alg_gauge; normalize=true)
+    Er = rightenv(envs, pos, ψ)
+    Hi = H[pos]
+    α = alg.noise
+    T = promote_type(scalartype(ψ), scalartype(Hi))
+    V = left_virtualspace(AC)
+    combiner, Vpert = _get_combiner(T, V, left_virtualspace(Hi))
+    
+    @plansor pert[l; r s] := α * (combiner')[l; li lh] * AC[li, si; ri] * Hi[lh, s; si, rh] * Er[ri rh; r]
+
+    AC = _transpose_tail(AC)
+    AC_expanded = catcodomain(AC, pert)
+
+    C, AR = right_gauge(AC_expanded, alg_gauge)
+    B = ψ.AL[pos-1]
+    AL = catdomain(B, zeros(T, codomain(B) ← Vpert))
+    
+    normalize && normalize!(C)
+    ψ.AC[pos] = (C, AR)
+    ψ.AC[pos - 1] = (AL, C)
+    return ψ
+end
\ No newline at end of file
diff --git a/src/algorithms/post_expand/post_expand.jl b/src/algorithms/post_expand/post_expand.jl
new file mode 100644
index 000000000..32fc6a34c
--- /dev/null
+++ b/src/algorithms/post_expand/post_expand.jl
@@ -0,0 +1,11 @@
+struct NoExpand <: Algorithm end
+
+_update_post_expand(::NoExpand, iter, ϵ) = NoExpand()
+
+function post_expand!(pos::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, ::NoExpand, envs, AC, alg_gauge; normalize::Bool = false)
+    return left_gauge!(ψ, pos, AC, alg_gauge; normalize)
+end
+
+function post_expand!(pos::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, ::NoExpand, envs, AC, alg_gauge; normalize::Bool = false)
+    return right_gauge!(ψ, pos, AC, alg_gauge; normalize)
+end

From 80ade6705705e480b504b50b4a7c112bcdafdcbc Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Mon, 20 Jul 2026 14:37:27 +0200
Subject: [PATCH 02/19] add tests for dmrg3s

---
 test/algorithms/groundstate.jl | 52 ++++++++++++++++++++++++++++++++++
 test/setup/testsetup.jl        | 31 ++++++++++++++++++++
 2 files changed, 83 insertions(+)

diff --git a/test/algorithms/groundstate.jl b/test/algorithms/groundstate.jl
index 605d45814..a496e57ba 100644
--- a/test/algorithms/groundstate.jl
+++ b/test/algorithms/groundstate.jl
@@ -115,6 +115,58 @@ verbosity_conv = 1
         @test dim(left_virtualspace(ψ, L ÷ 2)) == D
     end
 
+    @testset "DMRG3S" begin
+        # start from a small bond so the post-expansion is exercised, mirroring CBEDMRG above
+        Random.seed!(1234)
+        ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2))
+        v₀ = variance(ψ₀, H)
+        alg_post_expand = DMRG3S(0.1, ExponentialDecay(0.7))  # TODO: match final constructor API
+        trscheme = truncrank(D)
+
+        # test logging
+        ψ, envs, δ = find_groundstate(
+            ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_post_expand, trscheme)
+        )
+
+        ψ, envs, δ = find_groundstate(
+            ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 10, alg_post_expand, trscheme), envs
+        )
+        v = variance(ψ, H)
+
+        # test using low variance
+        @test sum(δ) ≈ 0 atol = 1.0e-3
+        @test v < v₀
+        @test v < 1.0e-2
+        # the bond should have grown to the truncation target
+        @test dim(left_virtualspace(ψ, L ÷ 2)) == D
+    end
+
+    @testset "DMRG3S escapes local minimum (Hubig et al. 2015, Sec. VII A)" begin
+        L_heis = 20
+        H_heis = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L=L_heis)
+
+        Random.seed!(1234)
+        ψ_bad = bad_initial_state(H_heis, L_heis)
+
+        ψ_stuck, envs_stuck, δ_stuck = find_groundstate(
+            ψ_bad, H_heis, DMRG(; verbosity = verbosity_conv, maxiter = 30)
+        )
+        E_stuck = real(expectation_value(ψ_stuck, H_heis, envs_stuck))
+
+        alg_post_expand = DMRG3S(0.1, ExponentialDecay(0.8))
+        ψ_escape, envs_escape, δ_escape = find_groundstate(
+            ψ_bad, H_heis, DMRG(;
+                verbosity = verbosity_conv, maxiter = 30,
+                alg_post_expand, trscheme = truncrank(20),
+            )
+        )
+        E_escape = real(expectation_value(ψ_escape, H_heis, envs_escape))
+
+        # paper reports E(α=0) = -6.35479, E(α≠0) = -8.6824724 for this L_heis = 20, S = 1/2 AFM setup
+        @test E_escape < E_stuck - 1.0
+        @test isapprox(E_escape, -8.6824724; atol = 1.0e-4)
+    end
+
     @testset "GradientGrassmann" begin
         ψ₀ = FiniteMPS(randn, ComplexF64, 10, ℙ^2, ℙ^D)
         v₀ = variance(ψ₀, H)
diff --git a/test/setup/testsetup.jl b/test/setup/testsetup.jl
index 1ea5e9f24..527f5212c 100644
--- a/test/setup/testsetup.jl
+++ b/test/setup/testsetup.jl
@@ -12,6 +12,7 @@ using LinearAlgebra: Diagonal
 using Combinatorics: permutations
 using TensorKitTensors.SpinOperators: S_x, S_y, S_z, S_x_S_x, S_y_S_y, S_z_S_z, S_exchange, S_plus_S_min, S_min_S_plus
 using TensorKitTensors.FermionOperators: f_plus_f_min, f_min_f_plus, f_plus_f_plus, f_min_f_min, f_num, f_hopping
+using Random: shuffle
 
 # exports
 export S_x, S_y, S_z
@@ -21,6 +22,7 @@ export force_planar
 export symm_mul_mpo
 export transverse_field_ising, heisenberg_XXX, bilinear_biquadratic_model, XY_model, kitaev_model
 export classical_ising_tensors, classical_ising, sixvertex
+export bad_initial_state
 
 # using TensorOperations
 
@@ -241,4 +243,33 @@ function sixvertex(; a = 1.0, b = 1.0, c = 1.0)
     return InfiniteMPO([permute(TensorMap(d, ℂ^2 ⊗ ℂ^2, ℂ^2 ⊗ ℂ^2), ((1, 2), (4, 3)))])
 end
 
+# Functions for testing DMRG3S
+function distinct_bitstrings(n::Int, n_up::Int, n_states::Int)
+    trivial_state = vcat(ones(Int, n_up), zeros(Int, n - n_up))
+    results = Set{Vector{Int}}()
+    while length(results) < n_states
+        push!(results, shuffle(trivial_state))
+    end
+    return collect(results)
+end
+
+function bad_initial_state(H, L; T = ComplexF64, n_states = 20, n_fixed = 3)
+    physd = U1Space(-1 // 2 => 1, 1 // 2 => 1)
+    n_free = L - n_fixed
+    n_up = (n_free - n_fixed) ÷ 2  # 7 of 17 up, so the fixed +3/2 tail nets to Sz_total = 0
+
+    configs = distinct_bitstrings(n_free, n_up, n_states)
+
+    return sum(configs) do bits
+        charges = [b == 1 ? 1 // 2 : -1 // 2 for b in bits]
+        append!(charges, fill(1 // 2, n_fixed))  # the 3 fixed, all-up sites
+
+        q = cumsum(charges)
+        @assert q[end] == 0 "configuration doesn't land in the Sz_total = 0 sector"
+        Vs = [U1Space(qi => 1) for qi in q[1:(end - 1)]]  # L-1 internal bonds only;
+
+        return normalize!(FiniteMPS(T, fill(physd, L), Vs))
+    end
+end
+
 end

From 212aaf554817a5881fdf21d4f75b696d9fcd288c Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Mon, 20 Jul 2026 14:38:14 +0200
Subject: [PATCH 03/19] Style: runic

---
 src/algorithms/post_expand/dmrg3s.jl | 18 +++++++++---------
 test/algorithms/groundstate.jl       |  2 +-
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/algorithms/post_expand/dmrg3s.jl b/src/algorithms/post_expand/dmrg3s.jl
index 9e613c814..16b87e8ef 100644
--- a/src/algorithms/post_expand/dmrg3s.jl
+++ b/src/algorithms/post_expand/dmrg3s.jl
@@ -26,20 +26,20 @@ function _get_combiner(::Type{T}, V1, V2) where {T}
     return isomorphism(T, (V1 ⊗ V2) ← Vf), Vf
 end
 
-function post_expand!(pos::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::DMRG3S, envs, AC, alg_gauge; normalize=true)
+function post_expand!(pos::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::DMRG3S, envs, AC, alg_gauge; normalize = true)
     El = leftenv(envs, pos, ψ)
     Hi = H[pos]
     α = alg.noise
     T = promote_type(scalartype(ψ), scalartype(Hi))
     V = right_virtualspace(AC)
     combiner, Vpert = _get_combiner(T, V, right_virtualspace(Hi))
-    
-    @plansor pert[-1 -2; -3] := α * El[-1 1; 2] * AC[2 3; 4] * Hi[1 -2; 3 5] * combiner[4, 5; -3]    
+
+    @plansor pert[-1 -2; -3] := α * El[-1 1; 2] * AC[2 3; 4] * Hi[1 -2; 3 5] * combiner[4, 5; -3]
 
     AC_expanded = catdomain(AC, pert)
 
     AL, C = left_gauge(AC_expanded, alg_gauge)
-    B = _transpose_tail(ψ.AR[pos+1])
+    B = _transpose_tail(ψ.AR[pos + 1])
     AR = _transpose_front(catcodomain(B, zeros(T, Vpert ← domain(B))))
 
     normalize && normalize!(C)
@@ -48,25 +48,25 @@ function post_expand!(pos::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::DM
     return ψ
 end
 
-function post_expand!(pos::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::DMRG3S, envs, AC, alg_gauge; normalize=true)
+function post_expand!(pos::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::DMRG3S, envs, AC, alg_gauge; normalize = true)
     Er = rightenv(envs, pos, ψ)
     Hi = H[pos]
     α = alg.noise
     T = promote_type(scalartype(ψ), scalartype(Hi))
     V = left_virtualspace(AC)
     combiner, Vpert = _get_combiner(T, V, left_virtualspace(Hi))
-    
+
     @plansor pert[l; r s] := α * (combiner')[l; li lh] * AC[li, si; ri] * Hi[lh, s; si, rh] * Er[ri rh; r]
 
     AC = _transpose_tail(AC)
     AC_expanded = catcodomain(AC, pert)
 
     C, AR = right_gauge(AC_expanded, alg_gauge)
-    B = ψ.AL[pos-1]
+    B = ψ.AL[pos - 1]
     AL = catdomain(B, zeros(T, codomain(B) ← Vpert))
-    
+
     normalize && normalize!(C)
     ψ.AC[pos] = (C, AR)
     ψ.AC[pos - 1] = (AL, C)
     return ψ
-end
\ No newline at end of file
+end
diff --git a/test/algorithms/groundstate.jl b/test/algorithms/groundstate.jl
index a496e57ba..6279551b0 100644
--- a/test/algorithms/groundstate.jl
+++ b/test/algorithms/groundstate.jl
@@ -143,7 +143,7 @@ verbosity_conv = 1
 
     @testset "DMRG3S escapes local minimum (Hubig et al. 2015, Sec. VII A)" begin
         L_heis = 20
-        H_heis = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L=L_heis)
+        H_heis = heisenberg_XXX(ComplexF64, U1Irrep; spin = 1 // 2, L = L_heis)
 
         Random.seed!(1234)
         ψ_bad = bad_initial_state(H_heis, L_heis)

From d1b22bd395d034f4654c6fd23155072382052961 Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Tue, 14 Jul 2026 17:31:20 +0200
Subject: [PATCH 04/19] feat: add NoiseSchedule composition and
 ExponentialDecay threshold
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Add `FunctionalSchedule`, wrapping an arbitrary `(noise, iter, ϵ) -> noise`
callable as a `NoiseSchedule`. Overload `∘` on `NoiseSchedule` to compose
two schedules into a `FunctionalSchedule`, following the usual function
composition convention (`s2` applied first, `s1` applied to its result).

Add a `threshold` keyword to `ExponentialDecay`, snapping the decayed noise
to exactly zero once it falls below it -- avoids running the (cheap but
non-free) expansion step indefinitely on a vanishingly small amplitude.
This also makes `ExponentialDecay ∘ Warmup` a natural way to combine a
hard iteration cutoff with smooth decay.
---
 src/MPSKit.jl                        |  2 +-
 src/algorithms/post_expand/dmrg3s.jl | 44 ++++++++++++++++++++++++++--
 2 files changed, 42 insertions(+), 4 deletions(-)

diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index 3240f9a25..bb7fef530 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -39,7 +39,7 @@ export TDVP, TDVP2, WI, WII, TaylorCluster
 export changebonds, changebonds!
 export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand
 export post_expand!, NoExpand
-export NoiseSchedule, ExponentialDecay, Warmup, DMRG3S
+export NoiseSchedule, FunctionalSchedule, ExponentialDecay, Warmup, DMRG3S
 export propagator
 export DynamicalDMRG, NaiveInvert, Jeckelmann
 export exact_diagonalization, fidelity_susceptibility
diff --git a/src/algorithms/post_expand/dmrg3s.jl b/src/algorithms/post_expand/dmrg3s.jl
index 16b87e8ef..86f02ff0c 100644
--- a/src/algorithms/post_expand/dmrg3s.jl
+++ b/src/algorithms/post_expand/dmrg3s.jl
@@ -1,14 +1,52 @@
 abstract type NoiseSchedule end
 
-struct ExponentialDecay <: NoiseSchedule
-    decay_rate::Float64
+"""
+    FunctionalSchedule(f) <: NoiseSchedule
+
+Wrap an arbitrary callable `f(noise, iter, ϵ) -> noise` as a [`NoiseSchedule`](@ref).
+Used internally by `∘` to compose schedules, but can also be constructed directly for
+ad hoc schedules that don't warrant their own named type.
+"""
+struct FunctionalSchedule{F} <: NoiseSchedule
+    f::F
+end
+(s::FunctionalSchedule)(noise, iter, ϵ) = s.f(noise, iter, ϵ)
+
+"""
+    s1::NoiseSchedule ∘ s2::NoiseSchedule
+
+Compose two noise schedules: `s2` is applied first, and its result is fed through `s1`,
+i.e. `(s1 ∘ s2)(noise, iter, ϵ) == s1(s2(noise, iter, ϵ), iter, ϵ)` — the same convention
+as function composition in `Base`.
+"""
+Base.:∘(s1::NoiseSchedule, s2::NoiseSchedule) =
+    FunctionalSchedule((noise, iter, ϵ) -> s1(s2(noise, iter, ϵ), iter, ϵ))
+
+"""
+    ExponentialDecay(decay_rate; threshold = 0.0)
+
+Noise schedule that shrinks geometrically: `noise -> noise * decay_rate^iter`, snapped
+to exactly zero once it falls below `threshold`. Use `decay_rate < 1` for a standalone
+[`DMRG3S`](@ref) run that gradually turns enrichment off as the state converges; a
+nonzero `threshold` avoids running the (cheap, but non-free) expansion step
+indefinitely on a noise amplitude too small to matter.
+"""
+struct ExponentialDecay{T<:Real} <: NoiseSchedule
+    decay_rate::T
+    threshold::T
+end
+ExponentialDecay(decay_rate, threshold) = ExponentialDecay(promote(decay_rate, threshold)...)
+ExponentialDecay(decay_rate; threshold = zero(decay_rate)) = ExponentialDecay(decay_rate, threshold)
+
+function (s::ExponentialDecay)(noise, iter, ϵ)
+    decayed = noise * s.decay_rate^iter
+    return decayed < s.threshold ? zero(decayed) : decayed
 end
 
 struct Warmup <: NoiseSchedule
     iters::Int
 end
 
-(s::ExponentialDecay)(noise, iter, ϵ) = noise * s.decay_rate^iter
 (s::Warmup)(noise, iter, ϵ) = iter ≤ s.iters ? noise : zero(noise)
 
 struct DMRG3S{N, S <: NoiseSchedule} <: Algorithm

From 3eafe4ac94a440e43684eed7c3ba3036d1b7edde Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Mon, 20 Jul 2026 15:45:38 +0200
Subject: [PATCH 05/19] wip: fold DMRG3S/NoExpand into alg_gauge, single
 dispatch point

---
 src/MPSKit.jl                             |  2 +-
 src/algorithms/groundstate/dmrg.jl        | 15 ++++++++-------
 src/algorithms/post_expand/dmrg3s.jl      | 21 ++++++++++++---------
 src/algorithms/post_expand/post_expand.jl | 17 ++++++++---------
 4 files changed, 29 insertions(+), 26 deletions(-)

diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index bb7fef530..26916d7a8 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -38,7 +38,7 @@ export time_evolve, timestep, timestep!, make_time_mpo
 export TDVP, TDVP2, WI, WII, TaylorCluster
 export changebonds, changebonds!
 export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand
-export post_expand!, NoExpand
+export NoExpand
 export NoiseSchedule, FunctionalSchedule, ExponentialDecay, Warmup, DMRG3S
 export propagator
 export DynamicalDMRG, NaiveInvert, Jeckelmann
diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index de8ac4bd8..c1cc15399 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -71,7 +71,7 @@ function local_update!(
         site, direction,
         ψ, O, alg::DMRG, envs,
         ϵ_global, ϵ_trunc, decay_rate,
-        timeroutput
+        alg_gauge, timeroutput
     )
     ϵ_local = calc_galerkin(site, ψ, O, ψ, envs)
 
@@ -88,7 +88,7 @@ function local_update!(
     end
 
     # 3. gauge
-    ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge!(ψ, site, direction, AC′, alg.alg_gauge; normalize = true)
+    ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge!(site, direction, ψ, O, envs, AC′, alg_gauge; normalize = true)
 
     # 4. bookkeeping: measured contraction factor per matvec, kept a strict contraction in (0, 1)
     decay_rate = clamp((first(info.normres) / ϵ_local)^(1 / max(1, info.numops)), 1.0e-3, 0.999)
@@ -137,7 +137,7 @@ function DMRG2(;
     alg_eigsolve′ = alg_eigsolve isa NamedTuple ?
         Defaults.alg_eigsolve(; adaptive = true, alg_eigsolve...) : alg_eigsolve
     # two-site DMRG always truncates the enlarged bond back down, so the gauge is a truncated SVD
-    alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)
+    alg_gauge = NoExpand(MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme))
     return DMRG2(tol, maxiter, verbosity, alg_eigsolve′, alg_gauge, finalize)
 end
 
@@ -145,7 +145,7 @@ function local_update!(
         pos, direction,
         ψ, O, alg::DMRG2, envs,
         ϵ_global, ϵ_trunc, decay_rate,
-        timeroutput
+        alg_gauge, timeroutput
     )
     Heff = @timeit timeroutput "AC2_hamiltonian" AC2_hamiltonian(pos, ψ, O, ψ, envs)
 
@@ -164,7 +164,7 @@ function local_update!(
 
     # 2. gauge: truncated SVD split back into single-site tensors and install;
     #           the discarded weight is the truncation error
-    ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge2!(ψ, pos, direction, newA2center, alg.alg_gauge; normalize = true)
+    ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge2!(site, direction, ψ, O, envs, newA2center, alg_gauge; normalize = true)
 
     # 3. bookkeeping: measured contraction factor per matvec, kept a strict contraction in (0, 1)
     decay_rate = clamp((first(info.normres) / ϵ_local)^(1 / max(1, info.numops)), 1.0e-3, 0.999)
@@ -203,6 +203,7 @@ function find_groundstate!(
     LoggingExtras.withlevel(; alg.verbosity) do
         @infov 2 loginit!(log, ϵ_global, expectation_value(ψ, H, envs))
         for iter in 1:(alg.maxiter)
+            alg_gauge = _update_alg_gauge(alg.alg_gauge, iter, ϵ_global)
             @timeit timeroutput "sweep" begin
                 # left-to-right
                 for pos in fwd
@@ -211,7 +212,7 @@ function find_groundstate!(
                         pos, Val(:right),
                         ψ, H, alg, envs,
                         ϵ_global, ϵ_truncs[pos], decay_rates[pos],
-                        timeroutput
+                        alg_gauge, timeroutput
                     )
                     ϵ_global = maximum(ϵ_locals)
                 end
@@ -223,7 +224,7 @@ function find_groundstate!(
                         pos, Val(:left),
                         ψ, H, alg, envs,
                         ϵ_global, ϵ_truncs[pos], decay_rates[pos],
-                        timeroutput
+                        alg_gauge, timeroutput
                     )
                     ϵ_global = maximum(ϵ_locals)
                 end
diff --git a/src/algorithms/post_expand/dmrg3s.jl b/src/algorithms/post_expand/dmrg3s.jl
index 86f02ff0c..7631d27c7 100644
--- a/src/algorithms/post_expand/dmrg3s.jl
+++ b/src/algorithms/post_expand/dmrg3s.jl
@@ -49,14 +49,17 @@ end
 
 (s::Warmup)(noise, iter, ϵ) = iter ≤ s.iters ? noise : zero(noise)
 
-struct DMRG3S{N, S <: NoiseSchedule} <: Algorithm
+struct DMRG3S{N, S <: NoiseSchedule, A} <: Algorithm
     noise::N
     schedule::S
+    alg_gauge::A
 end
 
-function _update_post_expand(alg::DMRG3S, iter, ϵ)
+DMRG3S(noise, schedule) = DMRG3S(noise, schedule, nothing)
+
+function _update_alg_gauge(alg::DMRG3S, iter, ϵ)
     noise = alg.schedule(alg.noise, iter, ϵ)
-    return iszero(noise) ? NoExpand() : DMRG3S(noise, alg.schedule)
+    return iszero(noise) ? NoExpand(alg.alg_gauge) : DMRG3S(noise, alg.schedule, alg.alg_gauge)
 end
 
 function _get_combiner(::Type{T}, V1, V2) where {T}
@@ -64,7 +67,7 @@ function _get_combiner(::Type{T}, V1, V2) where {T}
     return isomorphism(T, (V1 ⊗ V2) ← Vf), Vf
 end
 
-function post_expand!(pos::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::DMRG3S, envs, AC, alg_gauge; normalize = true)
+function gauge!(pos::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, envs, AC, alg::DMRG3S; normalize = true)
     El = leftenv(envs, pos, ψ)
     Hi = H[pos]
     α = alg.noise
@@ -76,17 +79,17 @@ function post_expand!(pos::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::DM
 
     AC_expanded = catdomain(AC, pert)
 
-    AL, C = left_gauge(AC_expanded, alg_gauge)
+    AL, C, ϵ = left_gauge(AC_expanded, alg.alg_gauge)
     B = _transpose_tail(ψ.AR[pos + 1])
     AR = _transpose_front(catcodomain(B, zeros(T, Vpert ← domain(B))))
 
     normalize && normalize!(C)
     ψ.AC[pos] = (AL, C)
     ψ.AC[pos + 1] = (C, AR)
-    return ψ
+    return ψ, ϵ
 end
 
-function post_expand!(pos::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::DMRG3S, envs, AC, alg_gauge; normalize = true)
+function gauge!(pos::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, envs, AC, alg::DMRG3S; normalize = true)
     Er = rightenv(envs, pos, ψ)
     Hi = H[pos]
     α = alg.noise
@@ -99,12 +102,12 @@ function post_expand!(pos::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, alg::DMR
     AC = _transpose_tail(AC)
     AC_expanded = catcodomain(AC, pert)
 
-    C, AR = right_gauge(AC_expanded, alg_gauge)
+    C, AR, ϵ = right_gauge(AC_expanded, alg.alg_gauge)
     B = ψ.AL[pos - 1]
     AL = catdomain(B, zeros(T, codomain(B) ← Vpert))
 
     normalize && normalize!(C)
     ψ.AC[pos] = (C, AR)
     ψ.AC[pos - 1] = (AL, C)
-    return ψ
+    return ψ, ϵ
 end
diff --git a/src/algorithms/post_expand/post_expand.jl b/src/algorithms/post_expand/post_expand.jl
index 32fc6a34c..37f3cba99 100644
--- a/src/algorithms/post_expand/post_expand.jl
+++ b/src/algorithms/post_expand/post_expand.jl
@@ -1,11 +1,10 @@
-struct NoExpand <: Algorithm end
-
-_update_post_expand(::NoExpand, iter, ϵ) = NoExpand()
-
-function post_expand!(pos::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, ::NoExpand, envs, AC, alg_gauge; normalize::Bool = false)
-    return left_gauge!(ψ, pos, AC, alg_gauge; normalize)
+struct NoExpand{A} <: Algorithm
+    alg_gauge::A
 end
 
-function post_expand!(pos::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, ::NoExpand, envs, AC, alg_gauge; normalize::Bool = false)
-    return right_gauge!(ψ, pos, AC, alg_gauge; normalize)
-end
+_update_alg_gauge(alg::NoExpand, iter, ϵ) = alg
+
+gauge!(pos::Int, direction, ψ::AbstractFiniteMPS, H, envs, AC, alg::NoExpand; normalize::Bool = false) = 
+    gauge!(ψ, pos, direction, AC, alg.alg_gauge; normalize)
+gauge2!(pos::Int, direction, ψ::AbstractFiniteMPS, H, envs, AC2, alg::NoExpand; normalize::Bool = false) = 
+    gauge2!(ψ, pos, direction, AC2, alg.alg_gauge; normalize)
\ No newline at end of file

From e38a99bd983adfcf7c4455c8caceb4151612b7e6 Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Mon, 20 Jul 2026 15:54:58 +0200
Subject: [PATCH 06/19] runic styling

---
 src/algorithms/post_expand/dmrg3s.jl      | 2 +-
 src/algorithms/post_expand/post_expand.jl | 6 +++---
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/algorithms/post_expand/dmrg3s.jl b/src/algorithms/post_expand/dmrg3s.jl
index 7631d27c7..13f496339 100644
--- a/src/algorithms/post_expand/dmrg3s.jl
+++ b/src/algorithms/post_expand/dmrg3s.jl
@@ -31,7 +31,7 @@ to exactly zero once it falls below `threshold`. Use `decay_rate < 1` for a stan
 nonzero `threshold` avoids running the (cheap, but non-free) expansion step
 indefinitely on a noise amplitude too small to matter.
 """
-struct ExponentialDecay{T<:Real} <: NoiseSchedule
+struct ExponentialDecay{T <: Real} <: NoiseSchedule
     decay_rate::T
     threshold::T
 end
diff --git a/src/algorithms/post_expand/post_expand.jl b/src/algorithms/post_expand/post_expand.jl
index 37f3cba99..b11f8b0c8 100644
--- a/src/algorithms/post_expand/post_expand.jl
+++ b/src/algorithms/post_expand/post_expand.jl
@@ -4,7 +4,7 @@ end
 
 _update_alg_gauge(alg::NoExpand, iter, ϵ) = alg
 
-gauge!(pos::Int, direction, ψ::AbstractFiniteMPS, H, envs, AC, alg::NoExpand; normalize::Bool = false) = 
+gauge!(pos::Int, direction, ψ::AbstractFiniteMPS, H, envs, AC, alg::NoExpand; normalize::Bool = false) =
     gauge!(ψ, pos, direction, AC, alg.alg_gauge; normalize)
-gauge2!(pos::Int, direction, ψ::AbstractFiniteMPS, H, envs, AC2, alg::NoExpand; normalize::Bool = false) = 
-    gauge2!(ψ, pos, direction, AC2, alg.alg_gauge; normalize)
\ No newline at end of file
+gauge2!(pos::Int, direction, ψ::AbstractFiniteMPS, H, envs, AC2, alg::NoExpand; normalize::Bool = false) =
+    gauge2!(ψ, pos, direction, AC2, alg.alg_gauge; normalize)

From c8b7984f935e85f8e67037c61fc34dd249aed823 Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Mon, 20 Jul 2026 16:18:47 +0200
Subject: [PATCH 07/19] fixup! wip: fold DMRG3S/NoExpand into alg_gauge, single
 dispatch point

fix typo in variable name
---
 src/algorithms/groundstate/dmrg.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index c1cc15399..58f846fc0 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -164,7 +164,7 @@ function local_update!(
 
     # 2. gauge: truncated SVD split back into single-site tensors and install;
     #           the discarded weight is the truncation error
-    ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge2!(site, direction, ψ, O, envs, newA2center, alg_gauge; normalize = true)
+    ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge2!(pos, direction, ψ, O, envs, newA2center, alg_gauge; normalize = true)
 
     # 3. bookkeeping: measured contraction factor per matvec, kept a strict contraction in (0, 1)
     decay_rate = clamp((first(info.normres) / ϵ_local)^(1 / max(1, info.numops)), 1.0e-3, 0.999)

From ec493e4f8e239cf5f6755ecf542066ecf97160a7 Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Mon, 20 Jul 2026 19:31:03 +0200
Subject: [PATCH 08/19] update DMRG interface

---
 src/algorithms/groundstate/dmrg.jl        | 79 +++++++++++++++++------
 src/algorithms/post_expand/dmrg3s.jl      | 45 ++++++++++++-
 src/algorithms/post_expand/post_expand.jl | 13 ++++
 3 files changed, 118 insertions(+), 19 deletions(-)

diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index 58f846fc0..ec79c2d53 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -2,6 +2,16 @@
 # center-move as in textbook single-site DMRG)
 _truncates(::MatrixAlgebraKit.AbstractAlgorithm) = false
 _truncates(::MatrixAlgebraKit.TruncatedAlgorithm) = true
+_truncates(alg::NoExpand) = _truncates(alg.alg_gauge)
+_truncates(alg::DMRG3S) = _truncates(alg.alg_gauge)
+
+# a no-truncation `trscheme` selects a (bond-preserving) QR gauge, anything else a truncated SVD
+_build_inner_gauge(trscheme, alg_svd, alg_orth) =
+    trscheme isa MatrixAlgebraKit.NoTruncation ? alg_orth :
+        MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)
+
+_expands(::NoExpand) = false
+_expands(::DMRG3S) = true
 
 """
 $(TYPEDEF)
@@ -9,17 +19,37 @@ $(TYPEDEF)
 Density Matrix Renormalization Group algorithm for finding the dominant eigenvector.
 
 Each site update is, in order: (1) an optional bond expansion (`alg_expand`), (2) a single-site
-eigensolve, and (3) a gauge step (`alg_gauge`). With the defaults (`alg_expand = nothing` and a
-non-truncating QR `alg_gauge`) this is textbook single-site DMRG, which cannot change the bond
-dimension. Setting `alg_expand` to a bond-expansion algorithm (e.g. [`OptimalExpand`](@ref),
-[`RandExpand`](@ref), [`SketchedExpand`](@ref)) enriches the bond with directions orthogonal to
-the current state ahead of each eigensolve, recovering Controlled Bond Expansion (CBE) DMRG; a
-truncating `alg_gauge` is then desirable to cut the enlarged bond back down.
-
-The gauge algorithm is selected in the keyword constructor from the `trscheme` argument: when it
-is `notrunc()` the gauge is a QR decomposition (`alg_orth`, [`Householder`](@extref
-MatrixAlgebraKit.Householder) by default), otherwise it is a truncated SVD (`alg_svd` with the
-given `trscheme`).
+eigensolve, and (3) a gauge step (`alg_gauge`). With the defaults (`alg_expand = nothing` and
+`alg_gauge = NoExpand()`, a non-truncating QR gauge) this is textbook single-site DMRG, which
+cannot change the bond dimension. Setting `alg_expand` to a bond-expansion algorithm (e.g.
+[`OptimalExpand`](@ref), [`RandExpand`](@ref), [`SketchedExpand`](@ref)) enriches the bond with
+directions orthogonal to the current state ahead of each eigensolve, recovering Controlled Bond
+Expansion (CBE) DMRG. Setting `alg_gauge` to a bond-expanding gauge algorithm (e.g. [`DMRG3S`](@ref))
+instead enriches the bond as part of the gauge step, after the eigensolve. Either way, a
+truncating gauge (see below) is then desirable to cut the enlarged bond back down.
+
+# Choosing the gauge
+
+By default, `alg_gauge` is built for you from `trscheme`/`alg_svd`/`alg_orth`: `trscheme =
+notrunc()` (the default) gives a QR decomposition (`alg_orth`, [`Householder`](@extref
+MatrixAlgebraKit.Householder) by default), any other `trscheme` gives a truncated SVD (`alg_svd`
+with that `trscheme`).
+
+```julia
+DMRG()                            # QR gauge, no truncation
+DMRG(; trscheme = truncdim(50))   # truncated SVD gauge
+```
+
+To use a bond-expanding gauge such as [`DMRG3S`](@ref), pass it directly as `alg_gauge`; `trscheme`
+etc. are still routed through to build the *inner* gauge it wraps, exactly as above:
+
+```julia
+DMRG(; alg_gauge = DMRG3S(0.1, ExponentialDecay(0.7)), trscheme = truncdim(50))
+```
+
+If `alg_gauge` is instead given with its inner gauge already set (e.g. `DMRG3S(0.1, sched,
+some_gauge)`), `trscheme`/`alg_svd`/`alg_orth` must be left at their defaults — passing both is an
+error, since it leaves two conflicting sources for the same setting.
 
 ## Fields
 
@@ -44,24 +74,37 @@ struct DMRG{A, F, E, G} <: Algorithm
     "algorithm used to expand the bond ahead of each local update, or `nothing` for none"
     alg_expand::E
 
-    "factorization used for the post-update gauge: a QR algorithm (no truncation) or a truncated SVD"
+    "gauge algorithm applied after each local update: `NoExpand` for a plain gauge step (a QR
+    algorithm with no truncation, or a truncated SVD), or an algorithm that additionally expands
+    the bond beforehand (e.g. [`DMRG3S`](@ref))"
     alg_gauge::G
 end
 function DMRG(;
         tol = Defaults.tol, maxiter = Defaults.maxiter, alg_eigsolve = (;),
         verbosity = Defaults.verbosity, finalize = Defaults._finalize,
-        alg_expand = nothing, trscheme = notrunc(),
+        alg_expand = nothing, alg_gauge = NoExpand(), trscheme = nothing,
         alg_svd = Defaults.alg_svd(), alg_orth = Defaults.alg_orth()
     )
     # single-site DMRG defaults to the per-bond adaptive controller (`AdaptiveKrylov`); pass
     # `alg_eigsolve = (; adaptive = false, ...)` to opt out (the splat overrides the default).
     alg_eigsolve′ = alg_eigsolve isa NamedTuple ?
         Defaults.alg_eigsolve(; adaptive = true, alg_eigsolve...) : alg_eigsolve
-    # a no-truncation `trscheme` selects a (bond-preserving) QR gauge, anything else a truncated SVD
-    alg_gauge = trscheme isa MatrixAlgebraKit.NoTruncation ? alg_orth :
-        MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)
-    if !isnothing(alg_expand) && !_truncates(alg_gauge)
-        @warn "DMRG with `alg_expand` but no truncation (`trscheme = notrunc()`): the bond dimension will grow unboundedly each sweep."
+    
+    inner_gauge = alg_gauge.alg_gauge
+    if isnothing(inner_gauge)
+        trscheme = something(trscheme, notrunc()) # enforce trscheme default here
+        inner_gauge = _build_inner_gauge(trscheme, alg_svd, alg_orth)
+        alg_gauge = set_alg_gauge(alg_gauge, inner_gauge)
+    else
+        isnothing(trscheme) || throw(ArgumentError(
+            "`trscheme` was given together with an `alg_gauge` that already carries its own " *
+            "gauge algorithm (e.g. `DMRG3S(noise, schedule, some_gauge)`); set the truncation " *
+            "via one or the other, not both."
+        ))
+    end
+
+    if (!isnothing(alg_expand) || _expands(alg_gauge)) && !_truncates(alg_gauge)
+        @warn "DMRG with a bond-expanding `alg_expand` and/or `alg_gauge` but no truncation (`trscheme = notrunc()`): the bond dimension will grow unboundedly each sweep."
     end
     return DMRG(tol, maxiter, verbosity, alg_eigsolve′, finalize, alg_expand, alg_gauge)
 end
diff --git a/src/algorithms/post_expand/dmrg3s.jl b/src/algorithms/post_expand/dmrg3s.jl
index 13f496339..b2b7ed664 100644
--- a/src/algorithms/post_expand/dmrg3s.jl
+++ b/src/algorithms/post_expand/dmrg3s.jl
@@ -1,3 +1,12 @@
+"""
+    NoiseSchedule
+
+Supertype for controlling how a bond-expansion algorithm's noise/perturbation amplitude
+evolves across DMRG sweeps. A schedule `s` is called as `s(noise, iter, ϵ) -> noise`, given
+the current noise amplitude, the outer iteration count, and the current global convergence
+error, and returns the amplitude to use for the next iteration. Schedules can be composed
+with `∘`; see [`ExponentialDecay`](@ref), [`Warmup`](@ref), [`FunctionalSchedule`](@ref).
+"""
 abstract type NoiseSchedule end
 
 """
@@ -43,12 +52,44 @@ function (s::ExponentialDecay)(noise, iter, ϵ)
     return decayed < s.threshold ? zero(decayed) : decayed
 end
 
+"""
+    Warmup(iters::Int)
+
+Noise schedule that keeps the noise amplitude constant for the first `iters` outer
+iterations, then drops it to exactly zero. Useful composed with a decaying schedule (e.g.
+`ExponentialDecay(0.7) ∘ Warmup(5)`) to hold the perturbation at full strength for a few
+sweeps before starting to taper it off.
+"""
 struct Warmup <: NoiseSchedule
     iters::Int
 end
 
 (s::Warmup)(noise, iter, ϵ) = iter ≤ s.iters ? noise : zero(noise)
 
+"""
+    DMRG3S(noise, schedule::NoiseSchedule)
+
+Gauge algorithm wrapper that, at every site update, injects a Hamiltonian-derived
+perturbation of the just-optimized tensor before gauging — the "strictly single-site DMRG with
+subspace expansion" scheme of Hubig, McCulloch, Schollwöck & Wolf,
+[Phys. Rev. B 91, 155115 (2015)](https://doi.org/10.1103/PhysRevB.91.155115). This lets
+single-site DMRG introduce basis states/quantum-number sectors absent from the initial
+state, helping it escape local minima that plain single-site DMRG can get stuck in.
+
+`noise` is the initial perturbation amplitude; `schedule` (see [`ExponentialDecay`](@ref),
+[`Warmup`](@ref)) controls how it evolves across outer iterations, and once it decays to
+exactly zero the gauge step reverts to a plain [`NoExpand`](@ref) for the remainder of the
+run. As with [`NoExpand`](@ref), the actual factorization used to gauge the enriched tensor
+is filled in by `DMRG`'s constructor, not supplied here directly — see `DMRG`'s docstring
+for the calling convention:
+
+```julia
+DMRG(; alg_gauge = DMRG3S(0.1, ExponentialDecay(0.7)), trscheme = truncdim(50))
+```
+
+A truncating `trscheme` is strongly recommended alongside `DMRG3S`, to cut the perturbed
+bond back down each sweep — `DMRG`'s constructor warns if none is given.
+"""
 struct DMRG3S{N, S <: NoiseSchedule, A} <: Algorithm
     noise::N
     schedule::S
@@ -57,6 +98,8 @@ end
 
 DMRG3S(noise, schedule) = DMRG3S(noise, schedule, nothing)
 
+set_alg_gauge(alg::DMRG3S, inner_gauge) = DMRG3S(alg.noise, alg.schedule, inner_gauge)
+
 function _update_alg_gauge(alg::DMRG3S, iter, ϵ)
     noise = alg.schedule(alg.noise, iter, ϵ)
     return iszero(noise) ? NoExpand(alg.alg_gauge) : DMRG3S(noise, alg.schedule, alg.alg_gauge)
@@ -110,4 +153,4 @@ function gauge!(pos::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, envs, AC, alg:
     ψ.AC[pos] = (C, AR)
     ψ.AC[pos - 1] = (AL, C)
     return ψ, ϵ
-end
+end
\ No newline at end of file
diff --git a/src/algorithms/post_expand/post_expand.jl b/src/algorithms/post_expand/post_expand.jl
index b11f8b0c8..f77c7cfe1 100644
--- a/src/algorithms/post_expand/post_expand.jl
+++ b/src/algorithms/post_expand/post_expand.jl
@@ -1,9 +1,22 @@
+"""
+    NoExpand()
+
+Gauge algorithm wrapper for a plain, unperturbed gauge step: `ψ.AC[pos]` is simply gauged
+into `ψ` with no bond enrichment. This is `DMRG`'s default `alg_gauge`. The actual
+factorization used (QR or truncated SVD) is filled in by `DMRG`'s constructor from its
+`trscheme`/`alg_svd`/`alg_orth` keywords; `NoExpand()` on its own is a placeholder and not
+meant to be used outside of `DMRG(...)`.
+"""
 struct NoExpand{A} <: Algorithm
     alg_gauge::A
 end
 
+NoExpand() = NoExpand(nothing)
+
 _update_alg_gauge(alg::NoExpand, iter, ϵ) = alg
 
+set_alg_gauge(alg::NoExpand, inner_gauge) = NoExpand(inner_gauge)
+
 gauge!(pos::Int, direction, ψ::AbstractFiniteMPS, H, envs, AC, alg::NoExpand; normalize::Bool = false) =
     gauge!(ψ, pos, direction, AC, alg.alg_gauge; normalize)
 gauge2!(pos::Int, direction, ψ::AbstractFiniteMPS, H, envs, AC2, alg::NoExpand; normalize::Bool = false) =

From 6da8a6bbdc248b649b5a3f58f453020bc0621cc3 Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Mon, 20 Jul 2026 19:31:16 +0200
Subject: [PATCH 09/19] update test to new interface

---
 test/algorithms/groundstate.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/test/algorithms/groundstate.jl b/test/algorithms/groundstate.jl
index 6279551b0..02e644a92 100644
--- a/test/algorithms/groundstate.jl
+++ b/test/algorithms/groundstate.jl
@@ -120,16 +120,16 @@ verbosity_conv = 1
         Random.seed!(1234)
         ψ₀ = FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^(D ÷ 2))
         v₀ = variance(ψ₀, H)
-        alg_post_expand = DMRG3S(0.1, ExponentialDecay(0.7))  # TODO: match final constructor API
+        alg_gauge = DMRG3S(0.1, ExponentialDecay(0.7))  # TODO: match final constructor API
         trscheme = truncrank(D)
 
         # test logging
         ψ, envs, δ = find_groundstate(
-            ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_post_expand, trscheme)
+            ψ₀, H, DMRG(; verbosity = verbosity_full, maxiter = 2, alg_gauge, trscheme)
         )
 
         ψ, envs, δ = find_groundstate(
-            ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 10, alg_post_expand, trscheme), envs
+            ψ, H, DMRG(; verbosity = verbosity_conv, maxiter = 10, alg_gauge, trscheme), envs
         )
         v = variance(ψ, H)
 
@@ -153,11 +153,11 @@ verbosity_conv = 1
         )
         E_stuck = real(expectation_value(ψ_stuck, H_heis, envs_stuck))
 
-        alg_post_expand = DMRG3S(0.1, ExponentialDecay(0.8))
+        alg_gauge = DMRG3S(0.1, ExponentialDecay(0.8))
         ψ_escape, envs_escape, δ_escape = find_groundstate(
             ψ_bad, H_heis, DMRG(;
                 verbosity = verbosity_conv, maxiter = 30,
-                alg_post_expand, trscheme = truncrank(20),
+                alg_gauge, trscheme = truncrank(20),
             )
         )
         E_escape = real(expectation_value(ψ_escape, H_heis, envs_escape))

From 09d5732807a21ab84189c121f9fb21b29559512e Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Mon, 20 Jul 2026 21:59:33 +0200
Subject: [PATCH 10/19] fix: update svd_trunc! call to use inner_alg_gauge for
 DMRG2 and DMRG

---
 src/algorithms/approximate/fvomps.jl |  2 +-
 src/algorithms/groundstate/dmrg.jl   | 18 +++++++++++-------
 src/algorithms/post_expand/dmrg3s.jl |  2 +-
 3 files changed, 13 insertions(+), 9 deletions(-)

diff --git a/src/algorithms/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl
index d5ffc6cc4..aacdacf2b 100644
--- a/src/algorithms/approximate/fvomps.jl
+++ b/src/algorithms/approximate/fvomps.jl
@@ -8,7 +8,7 @@ function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environment
             ϵ = 0.0
             for pos in [1:(length(ψ) - 1); (length(ψ) - 2):-1:1]
                 AC2′ = AC2_projection(pos, ψ, Oϕ, envs)
-                al, c, ar, = svd_trunc!(AC2′, alg.alg_gauge)
+                al, c, ar, = svd_trunc!(AC2′, inner_alg_gauge(alg))
 
                 AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
                 ϵ = max(ϵ, norm(al * c * ar - AC2) / norm(AC2))
diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index ec79c2d53..418bd789f 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -8,7 +8,7 @@ _truncates(alg::DMRG3S) = _truncates(alg.alg_gauge)
 # a no-truncation `trscheme` selects a (bond-preserving) QR gauge, anything else a truncated SVD
 _build_inner_gauge(trscheme, alg_svd, alg_orth) =
     trscheme isa MatrixAlgebraKit.NoTruncation ? alg_orth :
-        MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)
+    MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)
 
 _expands(::NoExpand) = false
 _expands(::DMRG3S) = true
@@ -89,18 +89,20 @@ function DMRG(;
     # `alg_eigsolve = (; adaptive = false, ...)` to opt out (the splat overrides the default).
     alg_eigsolve′ = alg_eigsolve isa NamedTuple ?
         Defaults.alg_eigsolve(; adaptive = true, alg_eigsolve...) : alg_eigsolve
-    
+
     inner_gauge = alg_gauge.alg_gauge
     if isnothing(inner_gauge)
         trscheme = something(trscheme, notrunc()) # enforce trscheme default here
         inner_gauge = _build_inner_gauge(trscheme, alg_svd, alg_orth)
         alg_gauge = set_alg_gauge(alg_gauge, inner_gauge)
     else
-        isnothing(trscheme) || throw(ArgumentError(
-            "`trscheme` was given together with an `alg_gauge` that already carries its own " *
-            "gauge algorithm (e.g. `DMRG3S(noise, schedule, some_gauge)`); set the truncation " *
-            "via one or the other, not both."
-        ))
+        isnothing(trscheme) || throw(
+            ArgumentError(
+                "`trscheme` was given together with an `alg_gauge` that already carries its own " *
+                    "gauge algorithm (e.g. `DMRG3S(noise, schedule, some_gauge)`); set the truncation " *
+                    "via one or the other, not both."
+            )
+        )
     end
 
     if (!isnothing(alg_expand) || _expands(alg_gauge)) && !_truncates(alg_gauge)
@@ -227,6 +229,8 @@ _num_updates(::DMRG2, ψ) = length(ψ) - 1
 _sweep_ranges(::DMRG, ψ) = (1:(length(ψ) - 1), length(ψ):-1:2)
 _sweep_ranges(::DMRG2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 2):-1:1)
 
+inner_alg_gauge(alg::Union{DMRG, DMRG2}) = alg.alg_gauge.alg_gauge
+
 function find_groundstate!(
         ψ::AbstractFiniteMPS, H, alg::Union{DMRG, DMRG2}, envs = environments(ψ, H, ψ)
     )
diff --git a/src/algorithms/post_expand/dmrg3s.jl b/src/algorithms/post_expand/dmrg3s.jl
index b2b7ed664..04be23711 100644
--- a/src/algorithms/post_expand/dmrg3s.jl
+++ b/src/algorithms/post_expand/dmrg3s.jl
@@ -153,4 +153,4 @@ function gauge!(pos::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, envs, AC, alg:
     ψ.AC[pos] = (C, AR)
     ψ.AC[pos - 1] = (AL, C)
     return ψ, ϵ
-end
\ No newline at end of file
+end

From e88778e661856c0054d45d70ec0eedb1b9f13dce Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Tue, 21 Jul 2026 16:49:54 +0200
Subject: [PATCH 11/19] remove type restriction in `fuser`

---
 src/utility/utility.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/utility/utility.jl b/src/utility/utility.jl
index 83f67aadc..ac8c3b957 100644
--- a/src/utility/utility.jl
+++ b/src/utility/utility.jl
@@ -143,7 +143,7 @@ function check_length(a, b...)
     return L
 end
 
-function fuser(::Type{TorA}, V1::S, V2::S) where {TorA, S <: IndexSpace}
+function fuser(::Type{TorA}, V1, V2) where {TorA}
     return isomorphism(TorA, fuse(V1 ⊗ V2), V1 ⊗ V2)
 end
 

From bef50192db7fa71aa20507bcbde6e705fe930134 Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Tue, 21 Jul 2026 23:56:20 +0200
Subject: [PATCH 12/19] construct combiner with `fuser` and reorder `pert`
 construction to mimic `MPO_AC_Hamiltonian`

---
 src/algorithms/post_expand/dmrg3s.jl | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/src/algorithms/post_expand/dmrg3s.jl b/src/algorithms/post_expand/dmrg3s.jl
index 04be23711..193512043 100644
--- a/src/algorithms/post_expand/dmrg3s.jl
+++ b/src/algorithms/post_expand/dmrg3s.jl
@@ -105,20 +105,17 @@ function _update_alg_gauge(alg::DMRG3S, iter, ϵ)
     return iszero(noise) ? NoExpand(alg.alg_gauge) : DMRG3S(noise, alg.schedule, alg.alg_gauge)
 end
 
-function _get_combiner(::Type{T}, V1, V2) where {T}
-    Vf = fuse(V1 ⊗ V2)
-    return isomorphism(T, (V1 ⊗ V2) ← Vf), Vf
-end
-
 function gauge!(pos::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, envs, AC, alg::DMRG3S; normalize = true)
     El = leftenv(envs, pos, ψ)
     Hi = H[pos]
     α = alg.noise
     T = promote_type(scalartype(ψ), scalartype(Hi))
     V = right_virtualspace(AC)
-    combiner, Vpert = _get_combiner(T, V, right_virtualspace(Hi))
 
-    @plansor pert[-1 -2; -3] := α * El[-1 1; 2] * AC[2 3; 4] * Hi[1 -2; 3 5] * combiner[4, 5; -3]
+    combiner = fuser(T, V, right_virtualspace(Hi))'
+    Vpert = only(domain(combiner))
+
+    @plansor pert[-1 -2; -3] := α * El[-1 5; 4] * AC[4 2; 1] * Hi[5 -2; 2 3] * combiner[1 3; -3]
 
     AC_expanded = catdomain(AC, pert)
 
@@ -138,11 +135,15 @@ function gauge!(pos::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, envs, AC, alg:
     α = alg.noise
     T = promote_type(scalartype(ψ), scalartype(Hi))
     V = left_virtualspace(AC)
-    combiner, Vpert = _get_combiner(T, V, left_virtualspace(Hi))
 
-    @plansor pert[l; r s] := α * (combiner')[l; li lh] * AC[li, si; ri] * Hi[lh, s; si, rh] * Er[ri rh; r]
+    combiner = fuser(T, V, left_virtualspace(Hi))
+    Vpert = only(codomain(combiner))
+    combiner = _transpose_front(combiner)
 
+    @plansor pert[-1 -2; -3] := α * combiner[-1 5; 4] * AC[4 2; 1] * Hi[5 -2; 2 3] * Er[1 3; -3]
     AC = _transpose_tail(AC)
+    pert = _transpose_tail(pert)
+
     AC_expanded = catcodomain(AC, pert)
 
     C, AR, ϵ = right_gauge(AC_expanded, alg.alg_gauge)

From d4a0f064f72b6709dfe190648ebc3a3e8b4bd11f Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Thu, 23 Jul 2026 18:38:11 +0200
Subject: [PATCH 13/19] move update alg_gauge to local_update!

---
 src/algorithms/groundstate/dmrg.jl | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index 418bd789f..3a8e292f8 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -116,7 +116,7 @@ function local_update!(
         site, direction,
         ψ, O, alg::DMRG, envs,
         ϵ_global, ϵ_trunc, decay_rate,
-        alg_gauge, timeroutput
+        iter, timeroutput
     )
     ϵ_local = calc_galerkin(site, ψ, O, ψ, envs)
 
@@ -132,6 +132,8 @@ function local_update!(
         fixedpoint(H_effective, ac_old, :SR, alg_eigsolve)
     end
 
+    alg_gauge = _update_alg_gauge(alg.alg_gauge, iter, ϵ_global)
+
     # 3. gauge
     ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge!(site, direction, ψ, O, envs, AC′, alg_gauge; normalize = true)
 
@@ -190,7 +192,7 @@ function local_update!(
         pos, direction,
         ψ, O, alg::DMRG2, envs,
         ϵ_global, ϵ_trunc, decay_rate,
-        alg_gauge, timeroutput
+        iter, timeroutput
     )
     Heff = @timeit timeroutput "AC2_hamiltonian" AC2_hamiltonian(pos, ψ, O, ψ, envs)
 
@@ -207,6 +209,8 @@ function local_update!(
         (newA2center, info)
     end
 
+    alg_gauge = _update_alg_gauge(alg.alg_gauge, iter, ϵ_global)
+
     # 2. gauge: truncated SVD split back into single-site tensors and install;
     #           the discarded weight is the truncation error
     ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge2!(pos, direction, ψ, O, envs, newA2center, alg_gauge; normalize = true)
@@ -250,7 +254,6 @@ function find_groundstate!(
     LoggingExtras.withlevel(; alg.verbosity) do
         @infov 2 loginit!(log, ϵ_global, expectation_value(ψ, H, envs))
         for iter in 1:(alg.maxiter)
-            alg_gauge = _update_alg_gauge(alg.alg_gauge, iter, ϵ_global)
             @timeit timeroutput "sweep" begin
                 # left-to-right
                 for pos in fwd
@@ -259,7 +262,7 @@ function find_groundstate!(
                         pos, Val(:right),
                         ψ, H, alg, envs,
                         ϵ_global, ϵ_truncs[pos], decay_rates[pos],
-                        alg_gauge, timeroutput
+                        iter, timeroutput
                     )
                     ϵ_global = maximum(ϵ_locals)
                 end
@@ -271,7 +274,7 @@ function find_groundstate!(
                         pos, Val(:left),
                         ψ, H, alg, envs,
                         ϵ_global, ϵ_truncs[pos], decay_rates[pos],
-                        alg_gauge, timeroutput
+                        iter, timeroutput
                     )
                     ϵ_global = maximum(ϵ_locals)
                 end

From 196239d52ecbc88685757341f845ad4e249da07e Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Thu, 23 Jul 2026 22:55:34 +0200
Subject: [PATCH 14/19] update BlockTensorKit compatibility to version 0.3.15

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 2a2d97ebc..24aa30455 100644
--- a/Project.toml
+++ b/Project.toml
@@ -36,7 +36,7 @@ MPSKitAdaptExt = "Adapt"
 [compat]
 Accessors = "0.1"
 Adapt = "4"
-BlockTensorKit = "0.3.14"
+BlockTensorKit = "0.3.15"
 Compat = "3.47, 4.10"
 DocStringExtensions = "0.9.3"
 HalfIntegers = "1.6.0"

From 05711eee7d45aaa15094fe403e6bf6b3c9f9a38e Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Thu, 23 Jul 2026 22:55:46 +0200
Subject: [PATCH 15/19] add _get_combiner function and update gauge! to use
 MPO_AC_Hamiltonian for perturbation calculation

---
 src/algorithms/post_expand/dmrg3s.jl | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/src/algorithms/post_expand/dmrg3s.jl b/src/algorithms/post_expand/dmrg3s.jl
index 193512043..b041dc582 100644
--- a/src/algorithms/post_expand/dmrg3s.jl
+++ b/src/algorithms/post_expand/dmrg3s.jl
@@ -105,6 +105,12 @@ function _update_alg_gauge(alg::DMRG3S, iter, ϵ)
     return iszero(noise) ? NoExpand(alg.alg_gauge) : DMRG3S(noise, alg.schedule, alg.alg_gauge)
 end
 
+# similar to `fuser` but the contracted leg is flattened
+function _get_combiner(::Type{TorA}, V1, V2) where {TorA}
+    Vprod = V1 ⊗ V2
+    return isomorphism(TorA, oplus(fuse(Vprod)), Vprod)
+end
+
 function gauge!(pos::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, envs, AC, alg::DMRG3S; normalize = true)
     El = leftenv(envs, pos, ψ)
     Hi = H[pos]
@@ -112,10 +118,11 @@ function gauge!(pos::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, envs, AC, alg
     T = promote_type(scalartype(ψ), scalartype(Hi))
     V = right_virtualspace(AC)
 
-    combiner = fuser(T, V, right_virtualspace(Hi))'
+    combiner = _get_combiner(T, V, right_virtualspace(Hi))'
     Vpert = only(domain(combiner))
 
-    @plansor pert[-1 -2; -3] := α * El[-1 5; 4] * AC[4 2; 1] * Hi[5 -2; 2 3] * combiner[1 3; -3]
+    mpo_ac = MPO_AC_Hamiltonian(El, Hi, combiner)
+    pert = α * mpo_ac(AC)
 
     AC_expanded = catdomain(AC, pert)
 
@@ -136,11 +143,12 @@ function gauge!(pos::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, envs, AC, alg:
     T = promote_type(scalartype(ψ), scalartype(Hi))
     V = left_virtualspace(AC)
 
-    combiner = fuser(T, V, left_virtualspace(Hi))
+    combiner = _get_combiner(T, V, left_virtualspace(Hi))
     Vpert = only(codomain(combiner))
     combiner = _transpose_front(combiner)
 
-    @plansor pert[-1 -2; -3] := α * combiner[-1 5; 4] * AC[4 2; 1] * Hi[5 -2; 2 3] * Er[1 3; -3]
+    mpo_ac = MPO_AC_Hamiltonian(combiner, Hi, Er)
+    pert = α * mpo_ac(AC)
     AC = _transpose_tail(AC)
     pert = _transpose_tail(pert)
 

From 3fb1809d585cf289148919576ffaac78986fecbc Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Fri, 24 Jul 2026 14:47:57 +0200
Subject: [PATCH 16/19] change arguments order in `gauge[2]!` function

---
 src/algorithms/groundstate/dmrg.jl        | 4 ++--
 src/algorithms/post_expand/dmrg3s.jl      | 4 ++--
 src/algorithms/post_expand/post_expand.jl | 4 ++--
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index 3a8e292f8..292effe35 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -135,7 +135,7 @@ function local_update!(
     alg_gauge = _update_alg_gauge(alg.alg_gauge, iter, ϵ_global)
 
     # 3. gauge
-    ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge!(site, direction, ψ, O, envs, AC′, alg_gauge; normalize = true)
+    ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge!(ψ, site, direction, O, envs, AC′, alg_gauge; normalize = true)
 
     # 4. bookkeeping: measured contraction factor per matvec, kept a strict contraction in (0, 1)
     decay_rate = clamp((first(info.normres) / ϵ_local)^(1 / max(1, info.numops)), 1.0e-3, 0.999)
@@ -213,7 +213,7 @@ function local_update!(
 
     # 2. gauge: truncated SVD split back into single-site tensors and install;
     #           the discarded weight is the truncation error
-    ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge2!(pos, direction, ψ, O, envs, newA2center, alg_gauge; normalize = true)
+    ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge2!(ψ, pos, direction, O, envs, newA2center, alg_gauge; normalize = true)
 
     # 3. bookkeeping: measured contraction factor per matvec, kept a strict contraction in (0, 1)
     decay_rate = clamp((first(info.normres) / ϵ_local)^(1 / max(1, info.numops)), 1.0e-3, 0.999)
diff --git a/src/algorithms/post_expand/dmrg3s.jl b/src/algorithms/post_expand/dmrg3s.jl
index b041dc582..4c3efbece 100644
--- a/src/algorithms/post_expand/dmrg3s.jl
+++ b/src/algorithms/post_expand/dmrg3s.jl
@@ -111,7 +111,7 @@ function _get_combiner(::Type{TorA}, V1, V2) where {TorA}
     return isomorphism(TorA, oplus(fuse(Vprod)), Vprod)
 end
 
-function gauge!(pos::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, envs, AC, alg::DMRG3S; normalize = true)
+function gauge!(ψ::AbstractFiniteMPS, pos::Int, ::Val{:right}, H, envs, AC, alg::DMRG3S; normalize = true)
     El = leftenv(envs, pos, ψ)
     Hi = H[pos]
     α = alg.noise
@@ -136,7 +136,7 @@ function gauge!(pos::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, envs, AC, alg
     return ψ, ϵ
 end
 
-function gauge!(pos::Int, ::Val{:left}, ψ::AbstractFiniteMPS, H, envs, AC, alg::DMRG3S; normalize = true)
+function gauge!(ψ::AbstractFiniteMPS, pos::Int, ::Val{:left}, H, envs, AC, alg::DMRG3S; normalize = true)
     Er = rightenv(envs, pos, ψ)
     Hi = H[pos]
     α = alg.noise
diff --git a/src/algorithms/post_expand/post_expand.jl b/src/algorithms/post_expand/post_expand.jl
index f77c7cfe1..d3f99f247 100644
--- a/src/algorithms/post_expand/post_expand.jl
+++ b/src/algorithms/post_expand/post_expand.jl
@@ -17,7 +17,7 @@ _update_alg_gauge(alg::NoExpand, iter, ϵ) = alg
 
 set_alg_gauge(alg::NoExpand, inner_gauge) = NoExpand(inner_gauge)
 
-gauge!(pos::Int, direction, ψ::AbstractFiniteMPS, H, envs, AC, alg::NoExpand; normalize::Bool = false) =
+gauge!(ψ::AbstractFiniteMPS, pos::Int, direction, H, envs, AC, alg::NoExpand; normalize::Bool = false) =
     gauge!(ψ, pos, direction, AC, alg.alg_gauge; normalize)
-gauge2!(pos::Int, direction, ψ::AbstractFiniteMPS, H, envs, AC2, alg::NoExpand; normalize::Bool = false) =
+gauge2!(ψ::AbstractFiniteMPS, pos::Int, direction, H, envs, AC2, alg::NoExpand; normalize::Bool = false) =
     gauge2!(ψ, pos, direction, AC2, alg.alg_gauge; normalize)

From ae4d889269726f9f8a5da8354731c8ccd82e4fa5 Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Fri, 24 Jul 2026 14:58:17 +0200
Subject: [PATCH 17/19] add reference for Hubig et al. 2015 in DMRG3S
 documentation

---
 docs/src/assets/mpskit.bib           | 15 +++++++++++++++
 src/algorithms/post_expand/dmrg3s.jl |  7 +++++--
 2 files changed, 20 insertions(+), 2 deletions(-)

diff --git a/docs/src/assets/mpskit.bib b/docs/src/assets/mpskit.bib
index 296f4194b..0b6da1846 100644
--- a/docs/src/assets/mpskit.bib
+++ b/docs/src/assets/mpskit.bib
@@ -614,3 +614,18 @@ @article{zhang2023
   url = {https://link.aps.org/doi/10.1103/PhysRevLett.130.151602},
   abstract = {We show that the Klein bottle entropy [H.-H. Tu, Phys. Rev. Lett. 119, 261603 (2017)] for conformal field theories perturbed by a relevant operator is a universal function of the dimensionless coupling constant. The universal scaling of the Klein bottle entropy near criticality provides an efficient approach to extract the scaling dimension of lattice operators via data collapse. As paradigmatic examples, we validate the universal scaling of the Klein bottle entropy for Ising and {$\mathbb{Z}$}3 parafermion conformal field theories with various perturbations using numerical simulation with continuous matrix product operator approach.}
 }
+
+@article{Hubig2015,
+  title = {Strictly single-site DMRG algorithm with subspace expansion},
+  author = {Hubig, C. and McCulloch, I. P. and Schollw\"ock, U. and Wolf, F. A.},
+  journal = {Phys. Rev. B},
+  volume = {91},
+  issue = {15},
+  pages = {155115},
+  numpages = {10},
+  year = {2015},
+  month = {Apr},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevB.91.155115},
+  url = {https://link.aps.org/doi/10.1103/PhysRevB.91.155115}
+}
diff --git a/src/algorithms/post_expand/dmrg3s.jl b/src/algorithms/post_expand/dmrg3s.jl
index 4c3efbece..83daedfee 100644
--- a/src/algorithms/post_expand/dmrg3s.jl
+++ b/src/algorithms/post_expand/dmrg3s.jl
@@ -71,8 +71,7 @@ end
 
 Gauge algorithm wrapper that, at every site update, injects a Hamiltonian-derived
 perturbation of the just-optimized tensor before gauging — the "strictly single-site DMRG with
-subspace expansion" scheme of Hubig, McCulloch, Schollwöck & Wolf,
-[Phys. Rev. B 91, 155115 (2015)](https://doi.org/10.1103/PhysRevB.91.155115). This lets
+subspace expansion" scheme. This lets
 single-site DMRG introduce basis states/quantum-number sectors absent from the initial
 state, helping it escape local minima that plain single-site DMRG can get stuck in.
 
@@ -89,6 +88,10 @@ DMRG(; alg_gauge = DMRG3S(0.1, ExponentialDecay(0.7)), trscheme = truncdim(50))
 
 A truncating `trscheme` is strongly recommended alongside `DMRG3S`, to cut the perturbed
 bond back down each sweep — `DMRG`'s constructor warns if none is given.
+
+## References
+
+* [Hubig et al. Phys. Rev. B 91, 155115 (2015)](@cite Hubig2015)
 """
 struct DMRG3S{N, S <: NoiseSchedule, A} <: Algorithm
     noise::N

From d895a63a32c8a75355bdd734ea189efc26eae34f Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Fri, 24 Jul 2026 15:59:15 +0200
Subject: [PATCH 18/19] Remove NoExpand algorithm

---
 src/MPSKit.jl                             |  1 -
 src/algorithms/groundstate/dmrg.jl        | 27 +++++++++++-----------
 src/algorithms/post_expand/dmrg3s.jl      |  7 +++---
 src/algorithms/post_expand/post_expand.jl | 28 ++++++-----------------
 4 files changed, 24 insertions(+), 39 deletions(-)

diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index 26916d7a8..5b0b2ba32 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -38,7 +38,6 @@ export time_evolve, timestep, timestep!, make_time_mpo
 export TDVP, TDVP2, WI, WII, TaylorCluster
 export changebonds, changebonds!
 export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand
-export NoExpand
 export NoiseSchedule, FunctionalSchedule, ExponentialDecay, Warmup, DMRG3S
 export propagator
 export DynamicalDMRG, NaiveInvert, Jeckelmann
diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index 292effe35..94926a990 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -2,7 +2,6 @@
 # center-move as in textbook single-site DMRG)
 _truncates(::MatrixAlgebraKit.AbstractAlgorithm) = false
 _truncates(::MatrixAlgebraKit.TruncatedAlgorithm) = true
-_truncates(alg::NoExpand) = _truncates(alg.alg_gauge)
 _truncates(alg::DMRG3S) = _truncates(alg.alg_gauge)
 
 # a no-truncation `trscheme` selects a (bond-preserving) QR gauge, anything else a truncated SVD
@@ -10,7 +9,7 @@ _build_inner_gauge(trscheme, alg_svd, alg_orth) =
     trscheme isa MatrixAlgebraKit.NoTruncation ? alg_orth :
     MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)
 
-_expands(::NoExpand) = false
+_expands(alg) = false
 _expands(::DMRG3S) = true
 
 """
@@ -20,13 +19,14 @@ Density Matrix Renormalization Group algorithm for finding the dominant eigenvec
 
 Each site update is, in order: (1) an optional bond expansion (`alg_expand`), (2) a single-site
 eigensolve, and (3) a gauge step (`alg_gauge`). With the defaults (`alg_expand = nothing` and
-`alg_gauge = NoExpand()`, a non-truncating QR gauge) this is textbook single-site DMRG, which
-cannot change the bond dimension. Setting `alg_expand` to a bond-expansion algorithm (e.g.
-[`OptimalExpand`](@ref), [`RandExpand`](@ref), [`SketchedExpand`](@ref)) enriches the bond with
-directions orthogonal to the current state ahead of each eigensolve, recovering Controlled Bond
-Expansion (CBE) DMRG. Setting `alg_gauge` to a bond-expanding gauge algorithm (e.g. [`DMRG3S`](@ref))
-instead enriches the bond as part of the gauge step, after the eigensolve. Either way, a
-truncating gauge (see below) is then desirable to cut the enlarged bond back down.
+`alg_gauge = nothing`, a non-truncating QR gauge derived from `trscheme = notrunc()`) this is
+textbook single-site DMRG, which cannot change the bond dimension. Setting `alg_expand` to a
+bond-expansion algorithm (e.g. [`OptimalExpand`](@ref), [`RandExpand`](@ref), [`SketchedExpand`](@ref))
+enriches the bond with directions orthogonal to the current state ahead of each eigensolve,
+recovering Controlled Bond Expansion (CBE) DMRG. Setting `alg_gauge` to a bond-expanding gauge
+algorithm (e.g. [`DMRG3S`](@ref)) instead enriches the bond as part of the gauge step, after the
+eigensolve. Either way, a truncating gauge (see below) is then desirable to cut the enlarged
+bond back down.
 
 # Choosing the gauge
 
@@ -82,7 +82,7 @@ end
 function DMRG(;
         tol = Defaults.tol, maxiter = Defaults.maxiter, alg_eigsolve = (;),
         verbosity = Defaults.verbosity, finalize = Defaults._finalize,
-        alg_expand = nothing, alg_gauge = NoExpand(), trscheme = nothing,
+        alg_expand = nothing, alg_gauge = nothing, trscheme = nothing,
         alg_svd = Defaults.alg_svd(), alg_orth = Defaults.alg_orth()
     )
     # single-site DMRG defaults to the per-bond adaptive controller (`AdaptiveKrylov`); pass
@@ -90,8 +90,7 @@ function DMRG(;
     alg_eigsolve′ = alg_eigsolve isa NamedTuple ?
         Defaults.alg_eigsolve(; adaptive = true, alg_eigsolve...) : alg_eigsolve
 
-    inner_gauge = alg_gauge.alg_gauge
-    if isnothing(inner_gauge)
+    if isnothing(alg_gauge) || isnothing(alg_gauge.alg_gauge)
         trscheme = something(trscheme, notrunc()) # enforce trscheme default here
         inner_gauge = _build_inner_gauge(trscheme, alg_svd, alg_orth)
         alg_gauge = set_alg_gauge(alg_gauge, inner_gauge)
@@ -184,7 +183,7 @@ function DMRG2(;
     alg_eigsolve′ = alg_eigsolve isa NamedTuple ?
         Defaults.alg_eigsolve(; adaptive = true, alg_eigsolve...) : alg_eigsolve
     # two-site DMRG always truncates the enlarged bond back down, so the gauge is a truncated SVD
-    alg_gauge = NoExpand(MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme))
+    alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)
     return DMRG2(tol, maxiter, verbosity, alg_eigsolve′, alg_gauge, finalize)
 end
 
@@ -233,7 +232,7 @@ _num_updates(::DMRG2, ψ) = length(ψ) - 1
 _sweep_ranges(::DMRG, ψ) = (1:(length(ψ) - 1), length(ψ):-1:2)
 _sweep_ranges(::DMRG2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 2):-1:1)
 
-inner_alg_gauge(alg::Union{DMRG, DMRG2}) = alg.alg_gauge.alg_gauge
+inner_alg_gauge(alg::Union{DMRG, DMRG2}) = alg_gauge(alg.alg_gauge)
 
 function find_groundstate!(
         ψ::AbstractFiniteMPS, H, alg::Union{DMRG, DMRG2}, envs = environments(ψ, H, ψ)
diff --git a/src/algorithms/post_expand/dmrg3s.jl b/src/algorithms/post_expand/dmrg3s.jl
index 83daedfee..9daea57cb 100644
--- a/src/algorithms/post_expand/dmrg3s.jl
+++ b/src/algorithms/post_expand/dmrg3s.jl
@@ -77,8 +77,8 @@ state, helping it escape local minima that plain single-site DMRG can get stuck
 
 `noise` is the initial perturbation amplitude; `schedule` (see [`ExponentialDecay`](@ref),
 [`Warmup`](@ref)) controls how it evolves across outer iterations, and once it decays to
-exactly zero the gauge step reverts to a plain [`NoExpand`](@ref) for the remainder of the
-run. As with [`NoExpand`](@ref), the actual factorization used to gauge the enriched tensor
+exactly zero the gauge step reverts to a plain gauge shift for the remainder of the
+run. The actual factorization used to gauge the enriched tensor
 is filled in by `DMRG`'s constructor, not supplied here directly — see `DMRG`'s docstring
 for the calling convention:
 
@@ -102,10 +102,11 @@ end
 DMRG3S(noise, schedule) = DMRG3S(noise, schedule, nothing)
 
 set_alg_gauge(alg::DMRG3S, inner_gauge) = DMRG3S(alg.noise, alg.schedule, inner_gauge)
+alg_gauge(alg) = alg.alg_gauge
 
 function _update_alg_gauge(alg::DMRG3S, iter, ϵ)
     noise = alg.schedule(alg.noise, iter, ϵ)
-    return iszero(noise) ? NoExpand(alg.alg_gauge) : DMRG3S(noise, alg.schedule, alg.alg_gauge)
+    return iszero(noise) ? alg.alg_gauge : DMRG3S(noise, alg.schedule, alg.alg_gauge)
 end
 
 # similar to `fuser` but the contracted leg is flattened
diff --git a/src/algorithms/post_expand/post_expand.jl b/src/algorithms/post_expand/post_expand.jl
index d3f99f247..5f8579c37 100644
--- a/src/algorithms/post_expand/post_expand.jl
+++ b/src/algorithms/post_expand/post_expand.jl
@@ -1,23 +1,9 @@
-"""
-    NoExpand()
+_update_alg_gauge(alg, iter, ϵ) = alg
 
-Gauge algorithm wrapper for a plain, unperturbed gauge step: `ψ.AC[pos]` is simply gauged
-into `ψ` with no bond enrichment. This is `DMRG`'s default `alg_gauge`. The actual
-factorization used (QR or truncated SVD) is filled in by `DMRG`'s constructor from its
-`trscheme`/`alg_svd`/`alg_orth` keywords; `NoExpand()` on its own is a placeholder and not
-meant to be used outside of `DMRG(...)`.
-"""
-struct NoExpand{A} <: Algorithm
-    alg_gauge::A
-end
+set_alg_gauge(::Nothing, inner_gauge) = inner_gauge
+alg_gauge(alg) = alg
 
-NoExpand() = NoExpand(nothing)
-
-_update_alg_gauge(alg::NoExpand, iter, ϵ) = alg
-
-set_alg_gauge(alg::NoExpand, inner_gauge) = NoExpand(inner_gauge)
-
-gauge!(ψ::AbstractFiniteMPS, pos::Int, direction, H, envs, AC, alg::NoExpand; normalize::Bool = false) =
-    gauge!(ψ, pos, direction, AC, alg.alg_gauge; normalize)
-gauge2!(ψ::AbstractFiniteMPS, pos::Int, direction, H, envs, AC2, alg::NoExpand; normalize::Bool = false) =
-    gauge2!(ψ, pos, direction, AC2, alg.alg_gauge; normalize)
+gauge!(ψ::AbstractFiniteMPS, pos::Int, direction, H, envs, AC, alg_gauge; normalize::Bool = false) =
+    gauge!(ψ, pos, direction, AC, alg_gauge; normalize)
+gauge2!(ψ::AbstractFiniteMPS, pos::Int, direction, H, envs, AC2, alg_gauge; normalize::Bool = false) =
+    gauge2!(ψ, pos, direction, AC2, alg_gauge; normalize)

From 0a64391a92d681562499215d63a22d0cbbb0a788 Mon Sep 17 00:00:00 2001
From: VinceNeede <v.bisogno3@studenti.unipi.it>
Date: Fri, 24 Jul 2026 16:47:56 +0200
Subject: [PATCH 19/19] fixup! Remove NoExpand algorithm

---
 src/algorithms/post_expand/dmrg3s.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/algorithms/post_expand/dmrg3s.jl b/src/algorithms/post_expand/dmrg3s.jl
index 9daea57cb..a35a11d29 100644
--- a/src/algorithms/post_expand/dmrg3s.jl
+++ b/src/algorithms/post_expand/dmrg3s.jl
@@ -102,7 +102,7 @@ end
 DMRG3S(noise, schedule) = DMRG3S(noise, schedule, nothing)
 
 set_alg_gauge(alg::DMRG3S, inner_gauge) = DMRG3S(alg.noise, alg.schedule, inner_gauge)
-alg_gauge(alg) = alg.alg_gauge
+alg_gauge(alg::DMRG3S) = alg.alg_gauge
 
 function _update_alg_gauge(alg::DMRG3S, iter, ϵ)
     noise = alg.schedule(alg.noise, iter, ϵ)