From c6290cd8aa92d4393987f309031e921f95a891d0 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/15] 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 1e36effa6afaf5e49c1880f99f155e7b16663e47 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/15] 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 02cd32cd3d1fc847ce747a3a08cf38673e9ae8b5 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/15] 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 df140147208e5c43e4216ee8ba03f4950932bf09 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/15] 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 103b0af3cee30aa4e269636818172f8c7a40bd4b 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/15] 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 80f413d6a0c5f53569972f44f8fce19002eda192 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/15] 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 cab229655a23ece98bd653a73c57cb5a66141ea7 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/15] 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 cdb150dc700e84541067218076c233612cba4d6f 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/15] 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 f50f2eb7c1c0220a60fc6fb238e5f4b7bd99a769 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/15] 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 e9c62ffb2d43e9072fbe56d03c408f8b4e366df4 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/15] 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 25909dd2be563c75a70a24c9d995aebe6e596b50 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/15] 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 451c745ad189271f9584563a1677ef3f280ccddf 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/15] 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 df6dbb206a39c2e7921a5ed61b7f8da023db237c 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/15] 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 35c340825549ce56338c120c82049a53c1e5a596 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/15] 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 660b14e0fdca3afa8e82e7cc1a42af4d01f3efdc 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/15] 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)