From da6174e733b42bd5e0e179e74be74e09172adcbd Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Thu, 23 Oct 2025 22:53:43 +0800
Subject: [PATCH 01/36] Improve time evolution interface

---
 examples/heisenberg_su/main.jl                |   6 +-
 src/PEPSKit.jl                                |   4 +-
 src/algorithms/time_evolution/evoltools.jl    |  19 +++
 src/algorithms/time_evolution/simpleupdate.jl | 142 +++++++-----------
 .../time_evolution/simpleupdate3site.jl       |  66 ++++----
 src/algorithms/time_evolution/time_evolve.jl  |  25 +++
 test/bondenv/benv_fu.jl                       |   4 +-
 test/examples/heisenberg.jl                   |   6 +-
 test/examples/j1j2_finiteT.jl                 |  17 +--
 test/examples/tf_ising_finiteT.jl             |  20 +--
 test/timeevol/cluster_projectors.jl           |  12 +-
 test/timeevol/sitedep_truncation.jl           |  14 +-
 12 files changed, 169 insertions(+), 166 deletions(-)
 create mode 100644 src/algorithms/time_evolution/time_evolve.jl

diff --git a/examples/heisenberg_su/main.jl b/examples/heisenberg_su/main.jl
index 18d266bb..ad30075b 100644
--- a/examples/heisenberg_su/main.jl
+++ b/examples/heisenberg_su/main.jl
@@ -70,12 +70,12 @@ fix a truncation error (if that can be reached by remaining below `Dbond`):
 
 dts = [1.0e-2, 1.0e-3, 4.0e-4]
 tols = [1.0e-6, 1.0e-8, 1.0e-8]
-maxiter = 10000
+nstep = 10000
 trscheme_peps = truncerr(1.0e-10) & truncdim(Dbond)
 
 for (dt, tol) in zip(dts, tols)
-    alg = SimpleUpdate(dt, tol, maxiter, trscheme_peps)
-    global peps, wts, = simpleupdate(peps, H, alg, wts; bipartite = true)
+    alg = SimpleUpdate(; tol, trscheme = trscheme_peps, bipartite = true, check_interval = 500)
+    global peps, wts, = time_evolve(peps, H, dt, nstep, alg, wts)
 end
 
 md"""
diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 9e1c24b7..25216ebc 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -70,6 +70,7 @@ include("algorithms/truncation/fullenv_truncation.jl")
 include("algorithms/truncation/bond_truncation.jl")
 
 include("algorithms/time_evolution/evoltools.jl")
+include("algorithms/time_evolution/time_evolve.jl")
 include("algorithms/time_evolution/simpleupdate.jl")
 include("algorithms/time_evolution/simpleupdate3site.jl")
 
@@ -99,7 +100,8 @@ export fixedpoint
 
 export absorb_weight
 export ALSTruncation, FullEnvTruncation
-export su_iter, su3site_iter, simpleupdate, SimpleUpdate
+export SimpleUpdate
+export time_evolve
 
 export InfiniteSquareNetwork
 export InfinitePartitionFunction
diff --git a/src/algorithms/time_evolution/evoltools.jl b/src/algorithms/time_evolution/evoltools.jl
index b153a488..ba36750e 100644
--- a/src/algorithms/time_evolution/evoltools.jl
+++ b/src/algorithms/time_evolution/evoltools.jl
@@ -1,5 +1,24 @@
 const InfiniteState = Union{InfinitePEPS, InfinitePEPO}
 
+function _process_timeevol_args(
+        state::InfiniteState, dt::Float64, imaginary_time::Bool, tol::Float64
+    )
+    dt′ = (state isa InfinitePEPS) ? dt : (dt / 2)
+    if (state isa InfinitePEPO)
+        @assert size(state)[3] == 1
+    end
+    if !imaginary_time
+        @assert (state isa InfinitePEPS) "Real time evolution of InfinitePEPO (Heisenberg picture) is not implemented."
+        dt′ = 1.0im * dt′
+    end
+    @assert tol >= 0
+    check_convergence = (tol > 0)
+    if check_convergence
+        @assert ((state isa InfinitePEPS) && imaginary_time) "Convergence check can only be enabled for ground state search."
+    end
+    return dt′, check_convergence
+end
+
 function MPSKit.infinite_temperature_density_matrix(H::LocalOperator)
     T = scalartype(H)
     A = map(physicalspace(H)) do Vp
diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index 78c5cd7b..c69ef2cb 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -1,24 +1,23 @@
 """
 $(TYPEDEF)
 
-Algorithm struct for simple update (SU) of infinite PEPS with bond weights.
-Each SU run is converged when the singular value difference becomes smaller than `tol`.
+Algorithm struct for simple update (SU) of InfinitePEPS or InfinitePEPO.
 
 ## Fields
 
 $(TYPEDFIELDS)
 """
-struct SimpleUpdate
-    dt::Number
-    tol::Float64
-    maxiter::Int
+@kwdef struct SimpleUpdate
     trscheme::TruncationScheme
+    check_interval::Int = 500
+    # only applicable to ground state search
+    bipartite::Bool = false
+    tol::Float64 = 0.0
+    # only applicable to PEPO evolution
+    gate_bothsides::Bool = false # when false, purified approach is assumed
 end
-# TODO: add kwarg constructor and SU Defaults
-
-"""
-$(SIGNATURES)
 
+#=
 Simple update of the x-bond between `[r,c]` and `[r,c+1]`.
 
 ```
@@ -26,7 +25,7 @@ Simple update of the x-bond between `[r,c]` and `[r,c+1]`.
     -- T[r,c] -- T[r,c+1] --
         |           |
 ```
-"""
+=#
 function _su_xbond!(
         state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
         row::Int, col::Int, trscheme::TruncationScheme; gate_ax::Int = 1
@@ -56,7 +55,7 @@ function _su_xbond!(
     return ϵ
 end
 
-"""
+#=
 Simple update of the y-bond between `[r,c]` and `[r-1,c]`.
 ```
         |
@@ -65,7 +64,7 @@ Simple update of the y-bond between `[r,c]` and `[r-1,c]`.
     -- T[r,c] ---
         |
 ```
-"""
+=#
 function _su_ybond!(
         state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
         row::Int, col::Int, trscheme::TruncationScheme; gate_ax::Int = 1
@@ -95,37 +94,24 @@ function _su_ybond!(
     return ϵ
 end
 
-"""
-    su_iter(state::InfinitePEPS, gate::LocalOperator, alg::SimpleUpdate, env::SUWeight; bipartite::Bool = false)
-    su_iter(densitymatrix::InfinitePEPO, gate::LocalOperator, alg::SimpleUpdate, env::SUWeight; gate_bothsides::Bool = true)
-
-One round of simple update, which applies the nearest neighbor `gate` on an InfinitePEPS `state` or InfinitePEPO `densitymatrix`.
-When the input `state` has a unit cell size of (2, 2), one can set `bipartite = true` to enforce the bipartite structure. 
-"""
 function su_iter(
-        state::InfiniteState, gate::LocalOperator, alg::SimpleUpdate, env::SUWeight;
-        bipartite::Bool = false, gate_bothsides::Bool = true
+        state::InfiniteState, gate::LocalOperator, alg::SimpleUpdate, env::SUWeight
     )
     @assert size(gate.lattice) == size(state)[1:2]
-    if state isa InfinitePEPS
-        gate_bothsides = false
-    else
-        @assert size(state, 3) == 1
-    end
     Nr, Nc = size(state)[1:2]
-    bipartite && (@assert Nr == Nc == 2)
+    alg.bipartite && (@assert Nr == Nc == 2)
     (Nr >= 2 && Nc >= 2) || throw(
         ArgumentError("`state` unit cell size for simple update should be no smaller than (2, 2)."),
     )
     state2, env2 = deepcopy(state), deepcopy(env)
-    gate_axs = gate_bothsides ? (1:2) : (1:1)
+    gate_axs = alg.gate_bothsides ? (1:2) : (1:1)
     for r in 1:Nr, c in 1:Nc
         term = get_gateterm(gate, (CartesianIndex(r, c), CartesianIndex(r, c + 1)))
         trscheme = truncation_scheme(alg.trscheme, 1, r, c)
         for gate_ax in gate_axs
             _su_xbond!(state2, term, env2, r, c, trscheme; gate_ax)
         end
-        if bipartite
+        if alg.bipartite
             rp1, cp1 = _next(r, Nr), _next(c, Nc)
             state2.A[rp1, cp1] = deepcopy(state2.A[r, c])
             state2.A[rp1, c] = deepcopy(state2.A[r, cp1])
@@ -136,92 +122,70 @@ function su_iter(
         for gate_ax in gate_axs
             _su_ybond!(state2, term, env2, r, c, trscheme; gate_ax)
         end
-        if bipartite
+        if alg.bipartite
             rm1, cm1 = _prev(r, Nr), _prev(c, Nc)
             state2.A[rm1, cm1] = deepcopy(state2.A[r, c])
             state2.A[r, cm1] = deepcopy(state2.A[rm1, c])
             env2.data[2, rm1, cm1] = deepcopy(env2.data[2, r, c])
         end
     end
-    return state2, env2
+    wtdiff = compare_weights(env2, env)
+    return state2, env2, (; wtdiff)
 end
 
 """
 Perform simple update with Hamiltonian `ham` containing up to nearest neighbor interaction terms. 
 """
 function _simpleupdate2site(
-        state::InfiniteState, ham::LocalOperator, alg::SimpleUpdate, env::SUWeight;
-        bipartite::Bool = false, check_interval::Int = 500, gate_bothsides::Bool = true
+        state::InfiniteState, H::LocalOperator,
+        dt::Float64, nstep::Int, alg::SimpleUpdate, env::SUWeight;
+        imaginary_time::Bool
     )
     time_start = time()
-    # exponentiating the 2-site Hamiltonian gate
-    if state isa InfinitePEPS
-        gate_bothsides = false
-    end
-    dt = gate_bothsides ? (alg.dt / 2) : alg.dt
-    gate = get_expham(ham, dt)
-    wtdiff = 1.0
-    env0 = deepcopy(env)
-    for count in 1:(alg.maxiter)
+    dt′, check_convergence = _process_timeevol_args(state, dt, imaginary_time, alg.tol)
+    gate = get_expham(H, dt′)
+    info = nothing
+    for step in 1:nstep
         time0 = time()
-        state, env = su_iter(state, gate, alg, env; bipartite, gate_bothsides)
-        wtdiff = compare_weights(env, env0)
-        converge = (wtdiff < alg.tol)
-        cancel = (count == alg.maxiter)
-        env0 = deepcopy(env)
+        state, env, info = su_iter(state, gate, alg, env)
+        converge = (info.wtdiff < alg.tol)
         time1 = time()
-        if ((count == 1) || (count % check_interval == 0) || converge || cancel)
+        if (step % alg.check_interval == 0) || (step == nstep) || converge
             @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))"
-            label = (converge ? "conv" : (cancel ? "cancel" : "iter"))
-            message = @sprintf(
-                "SU %s %-7d:  dt = %.0e,  weight diff = %.3e,  time = %.3f sec\n",
-                label, count, alg.dt, wtdiff, time1 - ((converge || cancel) ? time_start : time0)
+            @info @sprintf(
+                "SU iter %-7d: dt = %.0e, |Δλ| = %.3e. Time = %.3f s/it",
+                step, dt, info.wtdiff, time1 - time0
             )
-            cancel ? (@warn message) : (@info message)
         end
-        converge && break
+        if check_convergence
+            converge && break
+            if step == nstep
+                @warn "SU ground state search has not converged."
+            end
+        end
     end
-    return state, env, wtdiff
+    time_end = time()
+    @info @sprintf("Simple update finished. Time elasped: %.2f s", time_end - time_start)
+    return state, env, info
 end
 
-"""
-    simpleupdate(
-        state::InfinitePEPS, ham::LocalOperator, alg::SimpleUpdate, env::SUWeight;
-        bipartite::Bool = false, force_3site::Bool = false, check_interval::Int = 500
-    )
-    simpleupdate(
-        densitymatrix::InfinitePEPO, ham::LocalOperator, alg::SimpleUpdate, env::SUWeight;
-        gate_bothsides::Bool = true, force_3site::Bool = false, check_interval::Int = 500
-    )
-
-Perform a simple update on an InfinitePEPS `state` or an InfinitePEPO `densitymatrix`
-using the Hamiltonian `ham`, which can contain up to next-nearest-neighbor interaction terms.
-
-## Keyword Arguments
-
-- `bipartite::Bool=false`: If `true`, enforces the bipartite structure of the PEPS. 
-- `gate_bothsides::Bool=true`: (Effective only for PEPO evolution) If true, apply the Trotter gate to both side of the PEPO. Otherwise, the gate is only applied to the physical codomain legs. 
-- `force_3site::Bool=false`: Forces the use of the 3-site simple update algorithm, even if the Hamiltonian contains only nearest-neighbor terms.
-- `check_interval::Int=500`: Specifies the number of evolution steps between printing progress information.
-
-## Notes
-
-- Setting `bipartite = true` is allowed only for PEPS evolution with up to next-nearest neighbor terms, and requires the input `peps` to have a unit cell size of (2, 2). 
-"""
-function simpleupdate(
-        state::InfiniteState, ham::LocalOperator, alg::SimpleUpdate, env::SUWeight;
-        bipartite::Bool = false, gate_bothsides::Bool = true,
-        force_3site::Bool = false, check_interval::Int = 500
+function MPSKit.time_evolve(
+        state::InfiniteState, H::LocalOperator,
+        dt::Float64, nstep::Int, alg::SimpleUpdate, env::SUWeight;
+        imaginary_time::Bool = true, force_3site::Bool = false
     )
     # determine if Hamiltonian contains nearest neighbor terms only
-    nnonly = is_nearest_neighbour(ham)
+    nnonly = is_nearest_neighbour(H)
     use_3site = force_3site || !nnonly
-    @assert !(state isa InfinitePEPO && bipartite) "Evolution of PEPO with bipartite structure is not implemented."
-    @assert !(bipartite && use_3site) "3-site simple update is incompatible with bipartite lattice."
+    @assert !(alg.bipartite && state isa InfinitePEPO) "Evolution of PEPO with bipartite structure is not implemented."
+    @assert !(alg.bipartite && use_3site) "3-site simple update is incompatible with bipartite lattice."
+    if state isa InfinitePEPS
+        @assert !(alg.gate_bothsides) "alg.gate_bothsides = true is incompatible with PEPS evolution."
+    end
     # TODO: check SiteDependentTruncation is compatible with bipartite structure
     if use_3site
-        return _simpleupdate3site(state, ham, alg, env; gate_bothsides, check_interval)
+        return _simpleupdate3site(state, H, dt, nstep, alg, env; imaginary_time)
     else
-        return _simpleupdate2site(state, ham, alg, env; bipartite, gate_bothsides, check_interval)
+        return _simpleupdate2site(state, H, dt, nstep, alg, env; imaginary_time)
     end
 end
diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl
index 84cab37a..34c7062f 100644
--- a/src/algorithms/time_evolution/simpleupdate3site.jl
+++ b/src/algorithms/time_evolution/simpleupdate3site.jl
@@ -497,12 +497,9 @@ One round of 3-site simple update, which applies the Trotter gate MPOs `gatempos
 on an InfinitePEPS `state` or InfinitePEPO `densitymatrix`.
 """
 function su3site_iter(
-        state::InfiniteState, gatempos::Vector{G}, alg::SimpleUpdate, env::SUWeight;
-        gate_bothsides::Bool = true
+        state::InfiniteState, gatempos::Vector{G}, alg::SimpleUpdate, env::SUWeight
     ) where {G <: AbstractMatrix}
-    if state isa InfinitePEPS
-        gate_bothsides = false
-    else
+    if state isa InfinitePEPO
         @assert size(state, 3) == 1
     end
     Nr, Nc = size(state)[1:2]
@@ -521,53 +518,52 @@ function su3site_iter(
                 truncation_scheme(trscheme, 1, r, c)
                 truncation_scheme(trscheme, 2, r, _next(c, Nc))
             ]
-            _su3site_se!(state2, gs, env2, r, c, trschemes; gate_bothsides)
+            _su3site_se!(state2, gs, env2, r, c, trschemes; alg.gate_bothsides)
         end
         state2, env2 = rotl90(state2), rotl90(env2)
         trscheme = rotl90(trscheme)
     end
-    return state2, env2
+    wtdiff = compare_weights(env2, env)
+    return state2, env2, (; wtdiff)
 end
 
 """
-Perform 3-site simple update for Hamiltonian `ham`.
+Perform 3-site simple update for Hamiltonian `H`.
 """
 function _simpleupdate3site(
-        state::InfiniteState, ham::LocalOperator, alg::SimpleUpdate, env::SUWeight;
-        check_interval::Int = 500, gate_bothsides::Bool = true
+        state::InfiniteState, H::LocalOperator,
+        dt::Float64, nstep::Int, alg::SimpleUpdate, env::SUWeight;
+        imaginary_time::Bool
     )
     time_start = time()
-    # convert Hamiltonian to 3-site exponentiated gate MPOs
-    if state isa InfinitePEPS
-        gate_bothsides = false
-    end
-    dt = gate_bothsides ? (alg.dt / 2) : alg.dt
+    dt′, check_convergence = _process_timeevol_args(state, dt, imaginary_time, alg.tol)
     gatempos = [
-        _get_gatempos_se(ham, dt),
-        _get_gatempos_se(rotl90(ham), dt),
-        _get_gatempos_se(rot180(ham), dt),
-        _get_gatempos_se(rotr90(ham), dt),
+        _get_gatempos_se(H, dt′),
+        _get_gatempos_se(rotl90(H), dt′),
+        _get_gatempos_se(rot180(H), dt′),
+        _get_gatempos_se(rotr90(H), dt′),
     ]
-    wtdiff = 1.0
-    env0 = deepcopy(env)
-    for count in 1:(alg.maxiter)
+    info = nothing
+    for step in 1:nstep
         time0 = time()
-        state, env = su3site_iter(state, gatempos, alg, env; gate_bothsides)
-        wtdiff = compare_weights(env, env0)
-        converge = (wtdiff < alg.tol)
-        cancel = (count == alg.maxiter)
-        env0 = deepcopy(env)
+        state, env, info = su3site_iter(state, gatempos, alg, env)
+        converge = (info.wtdiff < alg.tol)
         time1 = time()
-        if ((count == 1) || (count % check_interval == 0) || converge || cancel)
+        if (step % alg.check_interval == 0) || (step == nstep) || converge
             @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))"
-            label = (converge ? "conv" : (cancel ? "cancel" : "iter"))
-            message = @sprintf(
-                "SU %s %-7d:  dt = %.0e,  weight diff = %.3e,  time = %.3f sec\n",
-                label, count, alg.dt, wtdiff, time1 - ((converge || cancel) ? time_start : time0)
+            @info @sprintf(
+                "SU iter %-7d: dt = %.0e, |Δλ| = %.3e. Time = %.3f s/it",
+                step, dt, info.wtdiff, time1 - time0
             )
-            cancel ? (@warn message) : (@info message)
         end
-        converge && break
+        if check_convergence
+            converge && break
+            if step == nstep
+                @warn "SU ground state search has not converged."
+            end
+        end
     end
-    return state, env, wtdiff
+    time_end = time()
+    @info @sprintf("Simple update finished. Time elasped: %.2f s", time_end - time_start)
+    return state, env, info
 end
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
new file mode 100644
index 00000000..62becbfc
--- /dev/null
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -0,0 +1,25 @@
+@doc """
+    time_evolve(ψ₀, H, dt, nstep, alg, envs₀; kwargs...) -> (ψ, envs, info)
+
+Time-evolve the initial state `ψ₀` with Hamiltonian `H` over 
+the time span `dt * nstep` based on Trotter decomposition. 
+
+## Arguments
+
+- `ψ₀::Union{InfinitePEPS, InfinitePEPO}`: Initial state.
+- `H::LocalOperator`: Hamiltonian operator (time-independent).
+- `dt::Float64`: Trotter time evolution step.
+- `nstep::Int`: Number of Trotter steps to be taken, such that the evolved time span is `dt * nstep`.
+- `alg`: Time evolution algorithm.
+- `envs₀`: Environment of the initial state.
+
+## Keyword Arguments
+
+- `imaginary_time::Bool=false`: if true, the time evolution is done with an imaginary time step
+    instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful for using this
+    function to compute the ground state of a Hamiltonian, or to compute finite-temperature
+    properties of a system.
+- `tol::Float64 = 0.0`: Tolerance of the change in SUWeight (for simple update; default 1.0e-9) 
+    or energy (for full update; default 1.0e-8) to determine if the ground state search has converged.
+"""
+time_evolve
diff --git a/test/bondenv/benv_fu.jl b/test/bondenv/benv_fu.jl
index d6a374a8..52bfca5c 100644
--- a/test/bondenv/benv_fu.jl
+++ b/test/bondenv/benv_fu.jl
@@ -13,8 +13,8 @@ function get_hubbard_state(t::Float64 = 1.0, U::Float64 = 8.0)
     Vphy = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 2, (1, 1 // 2) => 1, (1, -1 // 2) => 1)
     peps = InfinitePEPS(rand, ComplexF64, Vphy, Vphy; unitcell = (Nr, Nc))
     wts = SUWeight(peps)
-    alg = SimpleUpdate(1.0e-2, 1.0e-8, 10000, truncerr(1.0e-10) & truncdim(4))
-    peps, = simpleupdate(peps, H, alg, wts; bipartite = false, check_interval = 2000)
+    alg = SimpleUpdate(; tol = 1.0e-8, trscheme = truncerr(1.0e-10) & truncdim(4), check_interval = 2000)
+    peps, = time_evolve(peps, H, 1.0e-2, 10000, alg, wts)
     normalize!.(peps.A, Inf)
     return peps
 end
diff --git a/test/examples/heisenberg.jl b/test/examples/heisenberg.jl
index 9a390d91..ed8e1b11 100644
--- a/test/examples/heisenberg.jl
+++ b/test/examples/heisenberg.jl
@@ -66,12 +66,12 @@ end
     # simple update
     dts = [1.0e-2, 1.0e-3, 1.0e-3, 1.0e-4]
     tols = [1.0e-7, 1.0e-8, 1.0e-8, 1.0e-8]
-    maxiter = 5000
+    nstep = 5000
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
         Dbond2 = (n == 2) ? Dbond + 2 : Dbond
         trscheme = truncerr(1.0e-10) & truncdim(Dbond2)
-        alg = SimpleUpdate(dt, tol, maxiter, trscheme)
-        peps, wts, = simpleupdate(peps, ham, alg, wts; bipartite = false)
+        alg = SimpleUpdate(; trscheme, bipartite = false)
+        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts)
     end
 
     # measure physical quantities with CTMRG
diff --git a/test/examples/j1j2_finiteT.jl b/test/examples/j1j2_finiteT.jl
index b28ef7de..8b827a69 100644
--- a/test/examples/j1j2_finiteT.jl
+++ b/test/examples/j1j2_finiteT.jl
@@ -31,24 +31,23 @@ trscheme_pepo = truncdim(7) & truncerr(1.0e-12)
 check_interval = 100
 
 # PEPO approach
-dt, maxiter = 1.0e-3, 600
-alg = SimpleUpdate(dt, 0.0, maxiter, trscheme_pepo)
-pepo, wts, = simpleupdate(pepo0, ham, alg, wts0; check_interval, gate_bothsides = true)
+dt, nstep = 1.0e-3, 600
+alg = SimpleUpdate(; trscheme = trscheme_pepo, gate_bothsides = true, check_interval)
+pepo, wts, = time_evolve(pepo0, ham, dt, nstep, alg, wts0)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 energy = expectation_value(pepo, ham, env) / (Nr * Nc)
-@info "β = $(dt * maxiter): tr(ρH) = $(energy)"
+@info "β = $(dt * nstep): tr(ρH) = $(energy)"
 @test energy ≈ bm[2] atol = 5.0e-3
 
 # PEPS (purified PEPO) approach
-dt, maxiter = 1.0e-3, 300
-alg = SimpleUpdate(dt, 0.0, maxiter, trscheme_pepo)
-pepo, wts, = simpleupdate(pepo0, ham, alg, wts0; check_interval, gate_bothsides = false)
+alg = SimpleUpdate(; trscheme = trscheme_pepo, gate_bothsides = false, check_interval)
+pepo, wts, = time_evolve(pepo0, ham, dt, nstep, alg, wts0)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 energy = expectation_value(pepo, ham, env) / (Nr * Nc)
-@info "β = $(dt * maxiter): tr(ρH) = $(energy)"
+@info "β = $(dt * nstep) / 2: tr(ρH) = $(energy)"
 @test energy ≈ bm[1] atol = 5.0e-3
 
 env = converge_env(InfinitePEPS(pepo), 16)
 energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc)
-@info "β = 2 × $(dt * maxiter): ⟨ρ|H|ρ⟩ = $(energy)"
+@info "β = $(dt * nstep): ⟨ρ|H|ρ⟩ = $(energy)"
 @test energy ≈ bm[2] atol = 5.0e-3
diff --git a/test/examples/tf_ising_finiteT.jl b/test/examples/tf_ising_finiteT.jl
index d14f10ae..37b2484d 100644
--- a/test/examples/tf_ising_finiteT.jl
+++ b/test/examples/tf_ising_finiteT.jl
@@ -54,29 +54,29 @@ wts0 = SUWeight(pepo0)
 
 trscheme_pepo = truncdim(8) & truncerr(1.0e-12)
 
-dt, maxiter = 1.0e-3, 400
-β = dt * maxiter
-alg = SimpleUpdate(dt, 0.0, maxiter, trscheme_pepo)
+dt, nstep = 1.0e-3, 400
+β = dt * nstep
 
 # when g = 2, β = 0.4 and 2β = 0.8 belong to two phases (without and with nonzero σᶻ)
 
-# PEPO approach
-## results at β, or T = 2.5
-pepo, wts, = simpleupdate(pepo0, ham, alg, wts0; gate_bothsides = true)
+# PEPO approach: results at β, or T = 2.5
+alg = SimpleUpdate(; trscheme = trscheme_pepo, gate_bothsides = true)
+pepo, wts, = time_evolve(pepo0, ham, dt, nstep, alg, wts0)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 result_β = measure_mag(pepo, env)
 @info "Magnetization at T = $(1 / β)" result_β
 @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2)
 
-## results at 2β, or T = 1.25
-pepo, wts, = simpleupdate(pepo, ham, alg, wts; gate_bothsides = true)
+# continue to get results at 2β, or T = 1.25
+pepo, wts, = time_evolve(pepo, ham, dt, nstep, alg, wts)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 result_2β = measure_mag(pepo, env)
 @info "Magnetization at T = $(1 / (2β))" result_2β
 @test isapprox(abs.(result_2β), bm_2β, rtol = 1.0e-4)
 
-# purification approach (should match 2β result)
-pepo, = simpleupdate(pepo0, ham, alg, wts0; gate_bothsides = false)
+# Purification approach: results at 2β, or T = 1.25
+alg = SimpleUpdate(; trscheme = trscheme_pepo, gate_bothsides = false)
+pepo, = time_evolve(pepo0, ham, dt, 2 * nstep, alg, wts0)
 env = converge_env(InfinitePEPS(pepo), 8)
 result_2β′ = measure_mag(pepo, env; purified = true)
 @info "Magnetization at T = $(1 / (2β)) (purification approach)" result_2β′
diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl
index 2d3dc867..432d65c5 100644
--- a/test/timeevol/cluster_projectors.jl
+++ b/test/timeevol/cluster_projectors.jl
@@ -116,11 +116,11 @@ end
     # usual 2-site simple update, and measure energy
     dts = [1.0e-2, 1.0e-2, 5.0e-3]
     tols = [1.0e-8, 1.0e-8, 1.0e-8]
-    maxiter = 5000
+    nstep = 5000
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
         trscheme = truncerr(1.0e-10) & truncdim(n == 1 ? 4 : 2)
-        alg = SimpleUpdate(dt, tol, maxiter, trscheme)
-        peps, wts, = simpleupdate(peps, ham, alg, wts; bipartite = true, check_interval = 1000)
+        alg = SimpleUpdate(; tol, trscheme, bipartite = true, check_interval = 1000)
+        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts)
     end
     normalize!.(peps.A, Inf)
     env = CTMRGEnv(rand, Float64, peps, Espace)
@@ -132,10 +132,8 @@ end
     tols = [1.0e-8, 1.0e-8]
     trscheme = truncerr(1.0e-10) & truncdim(2)
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
-        alg = SimpleUpdate(dt, tol, maxiter, trscheme)
-        peps, wts, = simpleupdate(
-            peps, ham, alg, wts; check_interval = 1000, force_3site = true
-        )
+        alg = SimpleUpdate(; tol, trscheme, check_interval = 1000)
+        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts; force_3site = true)
     end
     normalize!.(peps.A, Inf)
     env, = leading_boundary(env, peps; tol = ctmrg_tol, trscheme = trscheme_env)
diff --git a/test/timeevol/sitedep_truncation.jl b/test/timeevol/sitedep_truncation.jl
index 13f3e1c1..3ee588d6 100644
--- a/test/timeevol/sitedep_truncation.jl
+++ b/test/timeevol/sitedep_truncation.jl
@@ -27,8 +27,8 @@ end
     # set trscheme to be compatible with bipartite structure
     bonddims = stack([[6 4; 4 6], [5 7; 7 5]]; dims = 1)
     trscheme = SiteDependentTruncation(collect(truncdim(d) for d in bonddims))
-    alg = SimpleUpdate(1.0e-2, 1.0e-14, 4, trscheme)
-    peps, env, = simpleupdate(peps0, ham, alg, env0; bipartite = true)
+    alg = SimpleUpdate(; trscheme, bipartite = true)
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
     # check bipartite structure is preserved
@@ -44,22 +44,22 @@ end
 
 @testset "Simple update: generic 2-site and 3-site" begin
     Nr, Nc = 3, 4
-    ham = real(heisenberg_XYZ(InfiniteSquare(Nr, Nc); Jx = 1.0, Jy = 1.0, Jz = 1.0))
     Random.seed!(100)
     peps0 = InfinitePEPS(rand, Float64, ℂ^2, ℂ^10; unitcell = (Nr, Nc))
     normalize!.(peps0.A, Inf)
     env0 = SUWeight(peps0)
     # Site dependent truncation
     bonddims = rand(2:8, 2, Nr, Nc)
-    @show bonddims
     trscheme = SiteDependentTruncation(collect(truncdim(d) for d in bonddims))
-    alg = SimpleUpdate(1.0e-2, 1.0e-14, 2, trscheme)
+    alg = SimpleUpdate(; trscheme)
     # 2-site SU
-    peps, env, = simpleupdate(peps0, ham, alg, env0; bipartite = false)
+    ham = real(heisenberg_XYZ(InfiniteSquare(Nr, Nc); Jx = 1.0, Jy = 1.0, Jz = 1.0))
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
     # 3-site SU
-    peps, env, = simpleupdate(peps0, ham, alg, env0; bipartite = false, force_3site = true)
+    ham = real(j1_j2_model(InfiniteSquare(Nr, Nc); J1 = 1.0, J2 = 0.5, sublattice = false))
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
 end

From 88ced8924e28685735eba05d75c8588cd7f7dd52 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Thu, 23 Oct 2025 22:57:14 +0800
Subject: [PATCH 02/36] Update docstring

---
 src/algorithms/time_evolution/time_evolve.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index 62becbfc..a88d87c5 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -19,7 +19,5 @@ the time span `dt * nstep` based on Trotter decomposition.
     instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful for using this
     function to compute the ground state of a Hamiltonian, or to compute finite-temperature
     properties of a system.
-- `tol::Float64 = 0.0`: Tolerance of the change in SUWeight (for simple update; default 1.0e-9) 
-    or energy (for full update; default 1.0e-8) to determine if the ground state search has converged.
 """
 time_evolve

From dd3d602b475efdfe704a2e8fb7eac546f09fcfeb Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Thu, 23 Oct 2025 23:33:58 +0800
Subject: [PATCH 03/36] Move `force_3site` into `SimpleUpdate`

---
 src/algorithms/time_evolution/simpleupdate.jl | 5 +++--
 test/timeevol/cluster_projectors.jl           | 4 ++--
 2 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index c69ef2cb..9c143719 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -10,6 +10,7 @@ $(TYPEDFIELDS)
 @kwdef struct SimpleUpdate
     trscheme::TruncationScheme
     check_interval::Int = 500
+    force_3site::Bool = false
     # only applicable to ground state search
     bipartite::Bool = false
     tol::Float64 = 0.0
@@ -172,11 +173,11 @@ end
 function MPSKit.time_evolve(
         state::InfiniteState, H::LocalOperator,
         dt::Float64, nstep::Int, alg::SimpleUpdate, env::SUWeight;
-        imaginary_time::Bool = true, force_3site::Bool = false
+        imaginary_time::Bool = true
     )
     # determine if Hamiltonian contains nearest neighbor terms only
     nnonly = is_nearest_neighbour(H)
-    use_3site = force_3site || !nnonly
+    use_3site = alg.force_3site || !nnonly
     @assert !(alg.bipartite && state isa InfinitePEPO) "Evolution of PEPO with bipartite structure is not implemented."
     @assert !(alg.bipartite && use_3site) "3-site simple update is incompatible with bipartite lattice."
     if state isa InfinitePEPS
diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl
index 432d65c5..90220aa8 100644
--- a/test/timeevol/cluster_projectors.jl
+++ b/test/timeevol/cluster_projectors.jl
@@ -132,8 +132,8 @@ end
     tols = [1.0e-8, 1.0e-8]
     trscheme = truncerr(1.0e-10) & truncdim(2)
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
-        alg = SimpleUpdate(; tol, trscheme, check_interval = 1000)
-        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts; force_3site = true)
+        alg = SimpleUpdate(; tol, trscheme, check_interval = 1000, force_3site = true)
+        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts)
     end
     normalize!.(peps.A, Inf)
     env, = leading_boundary(env, peps; tol = ctmrg_tol, trscheme = trscheme_env)

From 6e71ed928c518803be0552e31e01ce89939487d4 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 28 Oct 2025 21:06:17 +0800
Subject: [PATCH 04/36] Define `Base.iterate` for SimpleUpdate

---
 src/algorithms/time_evolution/evoltools.jl    |   9 +-
 src/algorithms/time_evolution/simpleupdate.jl | 157 ++++++++++--------
 .../time_evolution/simpleupdate3site.jl       | 124 +++++---------
 src/algorithms/time_evolution/time_evolve.jl  |  35 ++--
 4 files changed, 138 insertions(+), 187 deletions(-)

diff --git a/src/algorithms/time_evolution/evoltools.jl b/src/algorithms/time_evolution/evoltools.jl
index ba36750e..02c6f56a 100644
--- a/src/algorithms/time_evolution/evoltools.jl
+++ b/src/algorithms/time_evolution/evoltools.jl
@@ -1,7 +1,7 @@
 const InfiniteState = Union{InfinitePEPS, InfinitePEPO}
 
 function _process_timeevol_args(
-        state::InfiniteState, dt::Float64, imaginary_time::Bool, tol::Float64
+        state::InfiniteState, dt::Float64, imaginary_time::Bool
     )
     dt′ = (state isa InfinitePEPS) ? dt : (dt / 2)
     if (state isa InfinitePEPO)
@@ -11,12 +11,7 @@ function _process_timeevol_args(
         @assert (state isa InfinitePEPS) "Real time evolution of InfinitePEPO (Heisenberg picture) is not implemented."
         dt′ = 1.0im * dt′
     end
-    @assert tol >= 0
-    check_convergence = (tol > 0)
-    if check_convergence
-        @assert ((state isa InfinitePEPS) && imaginary_time) "Convergence check can only be enabled for ground state search."
-    end
-    return dt′, check_convergence
+    return dt′
 end
 
 function MPSKit.infinite_temperature_density_matrix(H::LocalOperator)
diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index 9c143719..cf19e616 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -7,14 +7,31 @@ Algorithm struct for simple update (SU) of InfinitePEPS or InfinitePEPO.
 
 $(TYPEDFIELDS)
 """
-@kwdef struct SimpleUpdate
+@kwdef struct SimpleUpdate <: TimeEvolution
+    # Initial state (InfinitePEPS or InfinitePEPO)
+    ψ0::InfiniteState
+    # Hamiltonian
+    H::LocalOperator
+    # Initial Bond weights
+    env0::SUWeight
+    # Trotter time step
+    dt::Float64
+    # number of iteration steps
+    nstep::Int
+    # Truncation scheme after applying Trotter gates
     trscheme::TruncationScheme
+    # Switch for imaginary or real time
+    imaginary_time::Bool = true
+    # controls interval to print information
     check_interval::Int = 500
+    # force usage of 3-site simple update
     force_3site::Bool = false
-    # only applicable to ground state search
+    # ---- only applicable to ground state search ----
+    # assume bipartite unit cell structure
     bipartite::Bool = false
+    # converged change of weight between two steps
     tol::Float64 = 0.0
-    # only applicable to PEPO evolution
+    # ---- only applicable to PEPO evolution ----
     gate_bothsides::Bool = false # when false, purified approach is assumed
 end
 
@@ -28,14 +45,14 @@ Simple update of the x-bond between `[r,c]` and `[r,c+1]`.
 ```
 =#
 function _su_xbond!(
-        state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
+        ψ::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
         row::Int, col::Int, trscheme::TruncationScheme; gate_ax::Int = 1
     ) where {T <: Number, S <: ElementarySpace}
-    Nr, Nc = size(state)[1:2]
+    Nr, Nc = size(ψ)[1:2]
     @assert 1 <= row <= Nr && 1 <= col <= Nc
     cp1 = _next(col, Nc)
     # absorb environment weights
-    A, B = state.A[row, col], state.A[row, cp1]
+    A, B = ψ.A[row, col], ψ.A[row, cp1]
     A = absorb_weight(A, env, row, col, (NORTH, SOUTH, WEST); inv = false)
     B = absorb_weight(B, env, row, cp1, (NORTH, SOUTH, EAST); inv = false)
     normalize!(A, Inf)
@@ -51,7 +68,7 @@ function _su_xbond!(
     normalize!(B, Inf)
     normalize!(s, Inf)
     # update tensor dict and weight on current bond
-    state.A[row, col], state.A[row, cp1] = A, B
+    ψ.A[row, col], ψ.A[row, cp1] = A, B
     env.data[1, row, col] = s
     return ϵ
 end
@@ -67,14 +84,14 @@ Simple update of the y-bond between `[r,c]` and `[r-1,c]`.
 ```
 =#
 function _su_ybond!(
-        state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
+        ψ::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
         row::Int, col::Int, trscheme::TruncationScheme; gate_ax::Int = 1
     ) where {T <: Number, S <: ElementarySpace}
-    Nr, Nc = size(state)[1:2]
+    Nr, Nc = size(ψ)[1:2]
     @assert 1 <= row <= Nr && 1 <= col <= Nc
     rm1 = _prev(row, Nr)
     # absorb environment weights
-    A, B = state.A[row, col], state.A[rm1, col]
+    A, B = ψ.A[row, col], ψ.A[rm1, col]
     A = absorb_weight(A, env, row, col, (EAST, SOUTH, WEST); inv = false)
     B = absorb_weight(B, env, rm1, col, (NORTH, EAST, WEST); inv = false)
     normalize!(A, Inf)
@@ -90,103 +107,101 @@ function _su_ybond!(
     normalize!(A, Inf)
     normalize!(B, Inf)
     normalize!(s, Inf)
-    state.A[row, col], state.A[rm1, col] = A, B
+    ψ.A[row, col], ψ.A[rm1, col] = A, B
     env.data[2, row, col] = s
     return ϵ
 end
 
 function su_iter(
-        state::InfiniteState, gate::LocalOperator, alg::SimpleUpdate, env::SUWeight
+        ψ::InfiniteState, gate::LocalOperator, alg::SimpleUpdate, env::SUWeight
     )
-    @assert size(gate.lattice) == size(state)[1:2]
-    Nr, Nc = size(state)[1:2]
+    @assert size(gate.lattice) == size(ψ)[1:2]
+    Nr, Nc = size(ψ)[1:2]
     alg.bipartite && (@assert Nr == Nc == 2)
     (Nr >= 2 && Nc >= 2) || throw(
         ArgumentError("`state` unit cell size for simple update should be no smaller than (2, 2)."),
     )
-    state2, env2 = deepcopy(state), deepcopy(env)
+    ψ2, env2 = deepcopy(ψ), deepcopy(env)
     gate_axs = alg.gate_bothsides ? (1:2) : (1:1)
     for r in 1:Nr, c in 1:Nc
         term = get_gateterm(gate, (CartesianIndex(r, c), CartesianIndex(r, c + 1)))
         trscheme = truncation_scheme(alg.trscheme, 1, r, c)
         for gate_ax in gate_axs
-            _su_xbond!(state2, term, env2, r, c, trscheme; gate_ax)
+            _su_xbond!(ψ2, term, env2, r, c, trscheme; gate_ax)
         end
         if alg.bipartite
             rp1, cp1 = _next(r, Nr), _next(c, Nc)
-            state2.A[rp1, cp1] = deepcopy(state2.A[r, c])
-            state2.A[rp1, c] = deepcopy(state2.A[r, cp1])
+            ψ2.A[rp1, cp1] = deepcopy(ψ2.A[r, c])
+            ψ2.A[rp1, c] = deepcopy(ψ2.A[r, cp1])
             env2.data[1, rp1, cp1] = deepcopy(env2.data[1, r, c])
         end
         term = get_gateterm(gate, (CartesianIndex(r, c), CartesianIndex(r - 1, c)))
         trscheme = truncation_scheme(alg.trscheme, 2, r, c)
         for gate_ax in gate_axs
-            _su_ybond!(state2, term, env2, r, c, trscheme; gate_ax)
+            _su_ybond!(ψ2, term, env2, r, c, trscheme; gate_ax)
         end
         if alg.bipartite
             rm1, cm1 = _prev(r, Nr), _prev(c, Nc)
-            state2.A[rm1, cm1] = deepcopy(state2.A[r, c])
-            state2.A[r, cm1] = deepcopy(state2.A[rm1, c])
+            ψ2.A[rm1, cm1] = deepcopy(ψ2.A[r, c])
+            ψ2.A[r, cm1] = deepcopy(ψ2.A[rm1, c])
             env2.data[2, rm1, cm1] = deepcopy(env2.data[2, r, c])
         end
     end
     wtdiff = compare_weights(env2, env)
-    return state2, env2, (; wtdiff)
+    return ψ2, env2, (; wtdiff)
 end
 
-"""
-Perform simple update with Hamiltonian `ham` containing up to nearest neighbor interaction terms. 
-"""
-function _simpleupdate2site(
-        state::InfiniteState, H::LocalOperator,
-        dt::Float64, nstep::Int, alg::SimpleUpdate, env::SUWeight;
-        imaginary_time::Bool
-    )
-    time_start = time()
-    dt′, check_convergence = _process_timeevol_args(state, dt, imaginary_time, alg.tol)
-    gate = get_expham(H, dt′)
-    info = nothing
-    for step in 1:nstep
-        time0 = time()
-        state, env, info = su_iter(state, gate, alg, env)
-        converge = (info.wtdiff < alg.tol)
-        time1 = time()
-        if (step % alg.check_interval == 0) || (step == nstep) || converge
-            @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))"
-            @info @sprintf(
-                "SU iter %-7d: dt = %.0e, |Δλ| = %.3e. Time = %.3f s/it",
-                step, dt, info.wtdiff, time1 - time0
-            )
-        end
-        if check_convergence
-            converge && break
-            if step == nstep
-                @warn "SU ground state search has not converged."
-            end
-        end
-    end
-    time_end = time()
-    @info @sprintf("Simple update finished. Time elasped: %.2f s", time_end - time_start)
-    return state, env, info
-end
-
-function MPSKit.time_evolve(
-        state::InfiniteState, H::LocalOperator,
-        dt::Float64, nstep::Int, alg::SimpleUpdate, env::SUWeight;
-        imaginary_time::Bool = true
-    )
-    # determine if Hamiltonian contains nearest neighbor terms only
-    nnonly = is_nearest_neighbour(H)
+# Initial iteration
+function Base.iterate(alg::SimpleUpdate)
+    # sanity checks
+    nnonly = is_nearest_neighbour(alg.H)
     use_3site = alg.force_3site || !nnonly
-    @assert !(alg.bipartite && state isa InfinitePEPO) "Evolution of PEPO with bipartite structure is not implemented."
+    @assert alg.tol >= 0
+    @assert !(alg.bipartite && alg.ψ0 isa InfinitePEPO) "Evolution of PEPO with bipartite structure is not implemented."
     @assert !(alg.bipartite && use_3site) "3-site simple update is incompatible with bipartite lattice."
-    if state isa InfinitePEPS
+    if alg.ψ0 isa InfinitePEPS
         @assert !(alg.gate_bothsides) "alg.gate_bothsides = true is incompatible with PEPS evolution."
     end
-    # TODO: check SiteDependentTruncation is compatible with bipartite structure
-    if use_3site
-        return _simpleupdate3site(state, H, dt, nstep, alg, env; imaginary_time)
+    # construct Trotter 2-site gates or 3-site gate-MPOs
+    dt′ = _process_timeevol_args(alg.ψ0, alg.dt, alg.imaginary_time)
+    gate = if use_3site
+        [
+            _get_gatempos_se(alg.H, dt′),
+            _get_gatempos_se(rotl90(alg.H), dt′),
+            _get_gatempos_se(rot180(alg.H), dt′),
+            _get_gatempos_se(rotr90(alg.H), dt′),
+        ]
     else
-        return _simpleupdate2site(state, H, dt, nstep, alg, env; imaginary_time)
+        get_expham(alg.H, dt′)
+    end
+    time0 = time()
+    ψ, env, info = su_iter(alg.ψ0, gate, alg, alg.env0)
+    elapsed_time = time() - time0
+    converged = false
+    return (ψ, env, info), (ψ, env, info, gate, converged, 1, elapsed_time)
+end
+
+# subsequent iterations
+function Base.iterate(alg::SimpleUpdate, state)
+    ψ0, env0, info, gate, converged, it, elapsed_time = state
+    check_convergence = (alg.tol > 0)
+    if (it % alg.check_interval == 0) || (it == 1) || (it == alg.nstep) || converged
+        @info "Space of x-weight at [1, 1] = $(space(env0[1, 1, 1], 1))"
+        @info @sprintf(
+            "SU iter %-7d: dt = %.0e, |Δλ| = %.3e. Time = %.3f s/it",
+            it, alg.dt, info.wtdiff, elapsed_time
+        )
+        if (it == alg.nstep)
+            check_convergence && (@warn "SU: bond weights have not converged.")
+            return nothing
+        elseif converged
+            @info "SU: bond weights have converged."
+            return nothing
+        end
     end
+    time0 = time()
+    ψ, env, info = su_iter(ψ0, gate, alg, env0)
+    elapsed_time = time() - time0
+    converged = check_convergence ? (info.wtdiff < alg.tol) : false
+    return (ψ, env, info), (ψ, env, info, gate, converged, it + 1, elapsed_time)
 end
diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl
index 34c7062f..abd37e4b 100644
--- a/src/algorithms/time_evolution/simpleupdate3site.jl
+++ b/src/algorithms/time_evolution/simpleupdate3site.jl
@@ -95,7 +95,7 @@ Here `-*-` is the twist on the physical axis.
 # Truncation of a bond on OBC-MPS
 
 Suppose we want to truncate the bond between 
-the n-th and the (n+1)-th sites such that the truncated state
+the n-th and the (n+1)-th sites such that the truncated ψ
 ```
     |ψ̃⟩  =  M[1]---...---M̃[n]---M̃[n+1]---...---M[N]
             ↓            ↓      ↓              ↓
@@ -125,14 +125,14 @@ Then the fidelity is just
     F(ψ̃) = (norm(s̃[n], 2) / norm(s[n], 2))^2
 ```
 =#
-"""
+#=
 Perform QR decomposition through a PEPS tensor
 ```
              ╱           ╱
     --R0----M---  →  ---Q--*-R1--
           ╱ |         ╱ |
 ```
-"""
+=#
 function qr_through(
         R0::MPSBondTensor, M::GenericMPSTensor{S, 4}; normalize::Bool = true
     ) where {S <: ElementarySpace}
@@ -151,14 +151,14 @@ function qr_through(
     return q, r
 end
 
-"""
+#=
 Perform LQ decomposition through a tensor
 ```
              ╱           ╱
     --L0-*--Q---  ←  ---M--*-L1--
           ╱ |         ╱ |
 ```
-"""
+=#
 function lq_through(
         M::GenericMPSTensor{S, 4}, L1::MPSBondTensor; normalize::Bool = true
     ) where {S <: ElementarySpace}
@@ -177,9 +177,9 @@ function lq_through(
     return l, q
 end
 
-"""
+#=
 Given a cluster `Ms`, find all `R`, `L` matrices on each internal bond
-"""
+=#
 function _get_allRLs(Ms::Vector{T}) where {T <: GenericMPSTensor{<:ElementarySpace, 4}}
     # M1 -- (R1,L1) -- M2 -- (R2,L2) -- M3
     N = length(Ms)
@@ -198,7 +198,7 @@ function _get_allRLs(Ms::Vector{T}) where {T <: GenericMPSTensor{<:ElementarySpa
     return Rs, Ls
 end
 
-"""
+#=
 Given the tensors `R`, `L` on a bond, construct 
 the projectors `Pa`, `Pb` and the new bond weight `s`
 such that the contraction of `Pa`, `s`, `Pb` is identity when `trunc = notrunc`,
@@ -211,7 +211,7 @@ The arrows between `Pa`, `s`, `Pb` are
     rev = true:  - Pa --→-- Pb - 
                     2 → s → 1
 ```
-"""
+=#
 function _proj_from_RL(
         r::MPSBondTensor, l::MPSBondTensor;
         trunc::TruncationScheme = notrunc(), rev::Bool = false,
@@ -227,10 +227,10 @@ function _proj_from_RL(
     return Pa, s, Pb, ϵ
 end
 
-"""
+#=
 Given a cluster `Ms` and the pre-calculated `R`, `L` bond matrices,
 find all projectors `Pa`, `Pb` and Schmidt weights `wts` on internal bonds.
-"""
+=#
 function _get_allprojs(
         Ms, Rs, Ls, trschemes::Vector{E}, revs::Vector{Bool}
     ) where {E <: TruncationScheme}
@@ -252,9 +252,9 @@ function _get_allprojs(
     return Pas, Pbs, wts, ϵs
 end
 
-"""
+#=
 Find projectors to truncate internal bonds of the cluster `Ms`.
-"""
+=#
 function _cluster_truncate!(
         Ms::Vector{T}, trschemes::Vector{E}, revs::Vector{Bool}
     ) where {T <: GenericMPSTensor{<:ElementarySpace, 4}, E <: TruncationScheme}
@@ -269,7 +269,7 @@ function _cluster_truncate!(
     return wts, ϵs, Pas, Pbs
 end
 
-"""
+#=
 Apply the gate MPO `gs` on the cluster `Ms`.
 When `gate_ax` is 1 or 2, the gate acts from the physical codomain or domain side.
 
@@ -293,7 +293,7 @@ In the cluster, the axes of each tensor use the MPS order
       4  2            5  2
     M[1 2 3 4; 5]  M[1 2 3 4 5; 6]
 ```
-"""
+=#
 function _apply_gatempo!(
         Ms::Vector{T1}, gs::Vector{T2}; gate_ax::Int = 1
     ) where {T1 <: GenericMPSTensor{<:ElementarySpace, 4}, T2 <: AbstractTensorMap}
@@ -423,7 +423,7 @@ const perms_se_pepo = map(invperms_se_pepo) do (p1, p2)
     p = invperm((p1..., p2...))
     return (p[1:(end - 1)], (p[end],))
 end
-"""
+#=
 Obtain the 3-site cluster in the "southeast corner" of a square plaquette.
 ``` 
     r-1         M3
@@ -432,29 +432,29 @@ Obtain the 3-site cluster in the "southeast corner" of a square plaquette.
     r   M1 -←- M2
         c      c+1
 ```
-"""
-function get_3site_se(state::InfiniteState, env::SUWeight, row::Int, col::Int)
-    Nr, Nc = size(state)
+=#
+function get_3site_se(ψ::InfiniteState, env::SUWeight, row::Int, col::Int)
+    Nr, Nc = size(ψ)
     rm1, cp1 = _prev(row, Nr), _next(col, Nc)
     coords_se = [(row, col), (row, cp1), (rm1, cp1)]
-    perms_se = isa(state, InfinitePEPS) ? perms_se_peps : perms_se_pepo
+    perms_se = isa(ψ, InfinitePEPS) ? perms_se_peps : perms_se_pepo
     Ms = map(zip(coords_se, perms_se, openaxs_se)) do (coord, perm, openaxs)
-        M = absorb_weight(state.A[CartesianIndex(coord)], env, coord[1], coord[2], openaxs)
+        M = absorb_weight(ψ.A[CartesianIndex(coord)], env, coord[1], coord[2], openaxs)
         return permute(M, perm)
     end
     return Ms
 end
 
 function _su3site_se!(
-        state::InfiniteState, gs::Vector{T}, env::SUWeight,
+        ψ::InfiniteState, gs::Vector{T}, env::SUWeight,
         row::Int, col::Int, trschemes::Vector{E};
         gate_bothsides::Bool = true
     ) where {T <: AbstractTensorMap, E <: TruncationScheme}
-    Nr, Nc = size(state)
+    Nr, Nc = size(ψ)
     @assert 1 <= row <= Nr && 1 <= col <= Nc
     rm1, cp1 = _prev(row, Nr), _next(col, Nc)
     # southwest 3-site cluster and arrow direction within it
-    Ms = get_3site_se(state, env, row, col)
+    Ms = get_3site_se(ψ, env, row, col)
     revs = [isdual(space(M, 1)) for M in Ms[2:end]]
     Vphys = [codomain(M, 2) for M in Ms]
     normalize!.(Ms, Inf)
@@ -467,103 +467,55 @@ function _su3site_se!(
     ϵs = nothing
     for gate_ax in gate_axs
         _apply_gatempo!(Ms, gs; gate_ax)
-        if isa(state, InfinitePEPO)
+        if isa(ψ, InfinitePEPO)
             Ms = [first(_fuse_physicalspaces(M)) for M in Ms]
         end
         wts, ϵs, = _cluster_truncate!(Ms, trschemes, revs)
-        if isa(state, InfinitePEPO)
+        if isa(ψ, InfinitePEPO)
             Ms = [first(_unfuse_physicalspace(M, Vphy)) for (M, Vphy) in zip(Ms, Vphys)]
         end
         for (wt, wt_idx) in zip(wts, wt_idxs)
             env[CartesianIndex(wt_idx)] = normalize(wt, Inf)
         end
     end
-    invperms_se = isa(state, InfinitePEPS) ? invperms_se_peps : invperms_se_pepo
+    invperms_se = isa(ψ, InfinitePEPS) ? invperms_se_peps : invperms_se_pepo
     for (M, coord, invperm, openaxs, Vphy) in zip(Ms, coords, invperms_se, openaxs_se, Vphys)
         # restore original axes order
         M = permute(M, invperm)
         # remove weights on open axes of the cluster
         M = absorb_weight(M, env, coord[1], coord[2], openaxs; inv = true)
-        state.A[CartesianIndex(coord)] = normalize(M, Inf)
+        ψ.A[CartesianIndex(coord)] = normalize(M, Inf)
     end
     return ϵs
 end
 
-"""
-    su3site_iter(state::InfinitePEPS, gatempos, alg::SimpleUpdate, env::SUWeight)
-    su3site_iter(densitymatrix::InfinitePEPO, gatempos, alg::SimpleUpdate, env::SUWeight; gate_bothsides::Bool = true)
-
-One round of 3-site simple update, which applies the Trotter gate MPOs `gatempos`
-on an InfinitePEPS `state` or InfinitePEPO `densitymatrix`.
-"""
-function su3site_iter(
-        state::InfiniteState, gatempos::Vector{G}, alg::SimpleUpdate, env::SUWeight
+function su_iter(
+        ψ::InfiniteState, gatempos::Vector{G}, alg::SimpleUpdate, env::SUWeight
     ) where {G <: AbstractMatrix}
-    if state isa InfinitePEPO
-        @assert size(state, 3) == 1
+    if ψ isa InfinitePEPO
+        @assert size(ψ, 3) == 1
     end
-    Nr, Nc = size(state)[1:2]
+    Nr, Nc = size(ψ)[1:2]
     (Nr >= 2 && Nc >= 2) || throw(
         ArgumentError(
             "iPEPS unit cell size for simple update should be no smaller than (2, 2)."
         ),
     )
-    state2, env2 = deepcopy(state), deepcopy(env)
+    ψ2, env2 = deepcopy(ψ), deepcopy(env)
     trscheme = alg.trscheme
     for i in 1:4
-        Nr, Nc = size(state2)[1:2]
+        Nr, Nc = size(ψ2)[1:2]
         for r in 1:Nr, c in 1:Nc
             gs = gatempos[i][r, c]
             trschemes = [
                 truncation_scheme(trscheme, 1, r, c)
                 truncation_scheme(trscheme, 2, r, _next(c, Nc))
             ]
-            _su3site_se!(state2, gs, env2, r, c, trschemes; alg.gate_bothsides)
+            _su3site_se!(ψ2, gs, env2, r, c, trschemes; alg.gate_bothsides)
         end
-        state2, env2 = rotl90(state2), rotl90(env2)
+        ψ2, env2 = rotl90(ψ2), rotl90(env2)
         trscheme = rotl90(trscheme)
     end
     wtdiff = compare_weights(env2, env)
-    return state2, env2, (; wtdiff)
-end
-
-"""
-Perform 3-site simple update for Hamiltonian `H`.
-"""
-function _simpleupdate3site(
-        state::InfiniteState, H::LocalOperator,
-        dt::Float64, nstep::Int, alg::SimpleUpdate, env::SUWeight;
-        imaginary_time::Bool
-    )
-    time_start = time()
-    dt′, check_convergence = _process_timeevol_args(state, dt, imaginary_time, alg.tol)
-    gatempos = [
-        _get_gatempos_se(H, dt′),
-        _get_gatempos_se(rotl90(H), dt′),
-        _get_gatempos_se(rot180(H), dt′),
-        _get_gatempos_se(rotr90(H), dt′),
-    ]
-    info = nothing
-    for step in 1:nstep
-        time0 = time()
-        state, env, info = su3site_iter(state, gatempos, alg, env)
-        converge = (info.wtdiff < alg.tol)
-        time1 = time()
-        if (step % alg.check_interval == 0) || (step == nstep) || converge
-            @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))"
-            @info @sprintf(
-                "SU iter %-7d: dt = %.0e, |Δλ| = %.3e. Time = %.3f s/it",
-                step, dt, info.wtdiff, time1 - time0
-            )
-        end
-        if check_convergence
-            converge && break
-            if step == nstep
-                @warn "SU ground state search has not converged."
-            end
-        end
-    end
-    time_end = time()
-    @info @sprintf("Simple update finished. Time elasped: %.2f s", time_end - time_start)
-    return state, env, info
+    return ψ2, env2, (; wtdiff)
 end
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index a88d87c5..7809e7fd 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -1,23 +1,12 @@
-@doc """
-    time_evolve(ψ₀, H, dt, nstep, alg, envs₀; kwargs...) -> (ψ, envs, info)
-
-Time-evolve the initial state `ψ₀` with Hamiltonian `H` over 
-the time span `dt * nstep` based on Trotter decomposition. 
-
-## Arguments
-
-- `ψ₀::Union{InfinitePEPS, InfinitePEPO}`: Initial state.
-- `H::LocalOperator`: Hamiltonian operator (time-independent).
-- `dt::Float64`: Trotter time evolution step.
-- `nstep::Int`: Number of Trotter steps to be taken, such that the evolved time span is `dt * nstep`.
-- `alg`: Time evolution algorithm.
-- `envs₀`: Environment of the initial state.
-
-## Keyword Arguments
-
-- `imaginary_time::Bool=false`: if true, the time evolution is done with an imaginary time step
-    instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful for using this
-    function to compute the ground state of a Hamiltonian, or to compute finite-temperature
-    properties of a system.
-"""
-time_evolve
+abstract type TimeEvolution end
+
+function MPSKit.time_evolve(alg::Alg) where {Alg <: TimeEvolution}
+    time_start = time()
+    result = nothing
+    for state in alg
+        result = state
+    end
+    time_end = time()
+    @info @sprintf("Simple update finished. Total time elasped: %.2f s", time_end - time_start)
+    return result
+end

From fe4cb7a7edf2407d7b86d5a64164084525f8c073 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 28 Oct 2025 21:28:48 +0800
Subject: [PATCH 05/36] Update tests and some examples

---
 examples/heisenberg_su/main.jl      |  7 +++-
 examples/j1j2_su/main.jl            | 16 +++++---
 test/examples/heisenberg.jl         | 16 ++++----
 test/examples/j1j2_finiteT.jl       | 34 +++++++++-------
 test/examples/tf_ising_finiteT.jl   | 60 +++++++++++++++++------------
 test/timeevol/cluster_projectors.jl | 26 ++++++++-----
 test/timeevol/sitedep_truncation.jl | 17 ++++----
 7 files changed, 103 insertions(+), 73 deletions(-)

diff --git a/examples/heisenberg_su/main.jl b/examples/heisenberg_su/main.jl
index ad30075b..2a917f0c 100644
--- a/examples/heisenberg_su/main.jl
+++ b/examples/heisenberg_su/main.jl
@@ -74,8 +74,11 @@ nstep = 10000
 trscheme_peps = truncerr(1.0e-10) & truncdim(Dbond)
 
 for (dt, tol) in zip(dts, tols)
-    alg = SimpleUpdate(; tol, trscheme = trscheme_peps, bipartite = true, check_interval = 500)
-    global peps, wts, = time_evolve(peps, H, dt, nstep, alg, wts)
+    alg = SimpleUpdate(;
+        ψ0 = peps, env0 = wts, H, dt, tol, nstep,
+        trscheme = trscheme_peps, bipartite = true, check_interval = 500
+    )
+    global peps, wts, = time_evolve(alg)
 end
 
 md"""
diff --git a/examples/j1j2_su/main.jl b/examples/j1j2_su/main.jl
index 8d1b0030..697d3aac 100644
--- a/examples/j1j2_su/main.jl
+++ b/examples/j1j2_su/main.jl
@@ -48,15 +48,17 @@ Therefore, we shall gradually increase $J_2 / J_1$ from 0.1 to 0.5, each time in
 on the previously evolved PEPS:
 """
 
-dt, tol, maxiter = 1.0e-2, 1.0e-8, 30000
+dt, tol, nstep = 1.0e-2, 1.0e-8, 30000
 check_interval = 4000
 trscheme_peps = truncerr(1.0e-10) & truncdim(Dbond)
-alg = SimpleUpdate(dt, tol, maxiter, trscheme_peps)
 for J2 in 0.1:0.1:0.5
     H = real( ## convert Hamiltonian `LocalOperator` to real floats
         j1_j2_model(ComplexF64, symm, InfiniteSquare(Nr, Nc); J1, J2, sublattice = false),
     )
-    global peps, wts, = simpleupdate(peps, H, alg, wts; check_interval)
+    alg = SimpleUpdate(;
+        ψ0 = peps, env0 = wts, H, dt, tol, nstep, trscheme = trscheme_peps, check_interval
+    )
+    global peps, wts, = time_evolve(alg)
 end
 
 md"""
@@ -69,8 +71,10 @@ tols = [1.0e-9, 1.0e-9]
 J2 = 0.5
 H = real(j1_j2_model(ComplexF64, symm, InfiniteSquare(Nr, Nc); J1, J2, sublattice = false))
 for (dt, tol) in zip(dts, tols)
-    alg′ = SimpleUpdate(dt, tol, maxiter, trscheme_peps)
-    global peps, wts, = simpleupdate(peps, H, alg′, wts; check_interval)
+    alg = SimpleUpdate(;
+        ψ0 = peps, env0 = wts, H, dt, tol, nstep, trscheme = trscheme_peps, check_interval
+    )
+    global peps, wts, = time_evolve(alg)
 end
 
 md"""
@@ -115,7 +119,7 @@ using MPSKit: randomize!
 noise_peps = InfinitePEPS(randomize!.(deepcopy(peps.A)))
 peps₀ = peps + 1.0e-1noise_peps
 peps_opt, env_opt, E_opt, = fixedpoint(
-    H, peps₀, env; optimizer_alg = (; tol = 1.0e-4, maxiter = 80)
+    H, peps₀, env; optimizer_alg = (; tol = 1.0e-4, nstep = 80)
 );
 
 md"""
diff --git a/test/examples/heisenberg.jl b/test/examples/heisenberg.jl
index 2dcec90d..37486053 100644
--- a/test/examples/heisenberg.jl
+++ b/test/examples/heisenberg.jl
@@ -59,10 +59,10 @@ end
     normalize!.(peps.A, Inf)
 
     # Heisenberg model Hamiltonian
-    ham = heisenberg_XYZ(InfiniteSquare(N1, N2); Jx = 1.0, Jy = 1.0, Jz = 1.0)
+    H = heisenberg_XYZ(InfiniteSquare(N1, N2); Jx = 1.0, Jy = 1.0, Jz = 1.0)
     # assert imaginary part is zero
-    @assert length(imag(ham).terms) == 0
-    ham = real(ham)
+    @assert length(imag(H).terms) == 0
+    H = real(H)
 
     # simple update
     dts = [1.0e-2, 1.0e-3, 1.0e-3, 1.0e-4]
@@ -71,23 +71,21 @@ end
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
         Dbond2 = (n == 2) ? Dbond + 2 : Dbond
         trscheme = truncerr(1.0e-10) & truncdim(Dbond2)
-        alg = SimpleUpdate(; trscheme, bipartite = false)
-        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts)
+        alg = SimpleUpdate(; ψ0 = peps, env0 = wts, H, dt, nstep, trscheme, bipartite = false)
+        peps, wts, = time_evolve(alg)
     end
 
     # measure physical quantities with CTMRG
     normalize!.(peps.A, Inf)
     env, = leading_boundary(CTMRGEnv(rand, Float64, peps, Espace), peps; tol = ctmrg_tol, maxiter = ctmrg_maxiter)
-    e_site = cost_function(peps, env, ham) / (N1 * N2)
+    e_site = cost_function(peps, env, H) / (N1 * N2)
     @info "Simple update energy = $e_site"
     # benchmark data from Phys. Rev. B 94, 035133 (2016)
     @test isapprox(e_site, -0.6594; atol = 1.0e-3)
 
     # continue with auto differentiation
     peps_final, env_final, E_final, = fixedpoint(
-        ham,
-        peps,
-        env;
+        H, peps, env;
         optimizer_alg = (; tol = gradtol, maxiter = 25),
         boundary_alg = (; maxiter = ctmrg_maxiter),
         gradient_alg = (; alg = :linsolver, solver_alg = (; alg = :gmres)),
diff --git a/test/examples/j1j2_finiteT.jl b/test/examples/j1j2_finiteT.jl
index 8b827a69..a23b3c56 100644
--- a/test/examples/j1j2_finiteT.jl
+++ b/test/examples/j1j2_finiteT.jl
@@ -20,34 +20,40 @@ function converge_env(state, χ::Int)
 end
 
 Nr, Nc = 2, 2
-ham = j1_j2_model(
+H = j1_j2_model(
     Float64, SU2Irrep, InfiniteSquare(Nr, Nc);
     J1 = 1.0, J2 = 0.5, sublattice = false
 )
-pepo0 = PEPSKit.infinite_temperature_density_matrix(ham)
-wts0 = SUWeight(pepo0)
+ψ0 = PEPSKit.infinite_temperature_density_matrix(H)
+wts0 = SUWeight(ψ0)
 # 7 = 1 (spin-0) + 2 x 3 (spin-1)
-trscheme_pepo = truncdim(7) & truncerr(1.0e-12)
+trscheme_ψ = truncdim(7) & truncerr(1.0e-12)
 check_interval = 100
 
 # PEPO approach
 dt, nstep = 1.0e-3, 600
-alg = SimpleUpdate(; trscheme = trscheme_pepo, gate_bothsides = true, check_interval)
-pepo, wts, = time_evolve(pepo0, ham, dt, nstep, alg, wts0)
-env = converge_env(InfinitePartitionFunction(pepo), 16)
-energy = expectation_value(pepo, ham, env) / (Nr * Nc)
+alg = SimpleUpdate(;
+    ψ0, env0 = wts0, H, dt, nstep, trscheme = trscheme_ψ,
+    gate_bothsides = true, check_interval
+)
+ψ, wts, = time_evolve(alg)
+env = converge_env(InfinitePartitionFunction(ψ), 16)
+energy = expectation_value(ψ, H, env) / (Nr * Nc)
 @info "β = $(dt * nstep): tr(ρH) = $(energy)"
 @test energy ≈ bm[2] atol = 5.0e-3
 
 # PEPS (purified PEPO) approach
-alg = SimpleUpdate(; trscheme = trscheme_pepo, gate_bothsides = false, check_interval)
-pepo, wts, = time_evolve(pepo0, ham, dt, nstep, alg, wts0)
-env = converge_env(InfinitePartitionFunction(pepo), 16)
-energy = expectation_value(pepo, ham, env) / (Nr * Nc)
+alg = SimpleUpdate(;
+    ψ0, env0 = wts0, H, dt, nstep, trscheme = trscheme_ψ,
+    gate_bothsides = false, check_interval
+)
+ψ, wts, = time_evolve(alg)
+env = converge_env(InfinitePartitionFunction(ψ), 16)
+energy = expectation_value(ψ, H, env) / (Nr * Nc)
 @info "β = $(dt * nstep) / 2: tr(ρH) = $(energy)"
 @test energy ≈ bm[1] atol = 5.0e-3
 
-env = converge_env(InfinitePEPS(pepo), 16)
-energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc)
+env = converge_env(InfinitePEPS(ψ), 16)
+energy = expectation_value(ψ, H, ψ, env) / (Nr * Nc)
 @info "β = $(dt * nstep): ⟨ρ|H|ρ⟩ = $(energy)"
 @test energy ≈ bm[2] atol = 5.0e-3
diff --git a/test/examples/tf_ising_finiteT.jl b/test/examples/tf_ising_finiteT.jl
index 37b2484d..b85c3835 100644
--- a/test/examples/tf_ising_finiteT.jl
+++ b/test/examples/tf_ising_finiteT.jl
@@ -23,36 +23,36 @@ function tfising_model(T::Type{<:Number}, lattice::InfiniteSquare; J = 1.0, g =
     return PEPSKit.nearest_neighbour_hamiltonian(spaces, gate)
 end
 
-function converge_env(state, χ::Int)
+function converge_env(ψ, χ::Int)
     trscheme1 = truncdim(4) & truncerr(1.0e-12)
-    env0 = CTMRGEnv(randn, Float64, state, ℂ^4)
-    env, = leading_boundary(env0, state; alg = :sequential, trscheme = trscheme1, tol = 1.0e-10)
+    env0 = CTMRGEnv(randn, Float64, ψ, ℂ^4)
+    env, = leading_boundary(env0, ψ; alg = :sequential, trscheme = trscheme1, tol = 1.0e-10)
     trscheme2 = truncdim(χ) & truncerr(1.0e-12)
-    env, = leading_boundary(env, state; alg = :sequential, trscheme = trscheme2, tol = 1.0e-10)
+    env, = leading_boundary(env, ψ; alg = :sequential, trscheme = trscheme2, tol = 1.0e-10)
     return env
 end
 
-function measure_mag(pepo::InfinitePEPO, env::CTMRGEnv; purified::Bool = false)
+function measure_mag(ψ::InfinitePEPO, env::CTMRGEnv; purified::Bool = false)
     r, c = 1, 1
-    lattice = physicalspace(pepo)
+    lattice = physicalspace(ψ)
     Mx = LocalOperator(lattice, ((r, c),) => σˣ(Float64, Trivial))
     Mz = LocalOperator(lattice, ((r, c),) => σᶻ(Float64, Trivial))
     if purified
-        magx = expectation_value(pepo, Mx, pepo, env)
-        magz = expectation_value(pepo, Mz, pepo, env)
+        magx = expectation_value(ψ, Mx, ψ, env)
+        magz = expectation_value(ψ, Mz, ψ, env)
     else
-        magx = expectation_value(pepo, Mx, env)
-        magz = expectation_value(pepo, Mz, env)
+        magx = expectation_value(ψ, Mx, env)
+        magz = expectation_value(ψ, Mz, env)
     end
     return [magx, magz]
 end
 
 Nr, Nc = 2, 2
-ham = tfising_model(Float64, InfiniteSquare(Nr, Nc); J = 1.0, g = 2.0)
-pepo0 = PEPSKit.infinite_temperature_density_matrix(ham)
-wts0 = SUWeight(pepo0)
+H = tfising_model(Float64, InfiniteSquare(Nr, Nc); J = 1.0, g = 2.0)
+ψ0 = PEPSKit.infinite_temperature_density_matrix(H)
+wts0 = SUWeight(ψ0)
 
-trscheme_pepo = truncdim(8) & truncerr(1.0e-12)
+trscheme_ψ = truncdim(8) & truncerr(1.0e-12)
 
 dt, nstep = 1.0e-3, 400
 β = dt * nstep
@@ -60,24 +60,34 @@ dt, nstep = 1.0e-3, 400
 # when g = 2, β = 0.4 and 2β = 0.8 belong to two phases (without and with nonzero σᶻ)
 
 # PEPO approach: results at β, or T = 2.5
-alg = SimpleUpdate(; trscheme = trscheme_pepo, gate_bothsides = true)
-pepo, wts, = time_evolve(pepo0, ham, dt, nstep, alg, wts0)
-env = converge_env(InfinitePartitionFunction(pepo), 16)
-result_β = measure_mag(pepo, env)
+alg = SimpleUpdate(;
+    ψ0, env0 = wts0, H, dt, nstep,
+    trscheme = trscheme_ψ, gate_bothsides = true
+)
+ψ, wts, = time_evolve(alg)
+env = converge_env(InfinitePartitionFunction(ψ), 16)
+result_β = measure_mag(ψ, env)
 @info "Magnetization at T = $(1 / β)" result_β
 @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2)
 
 # continue to get results at 2β, or T = 1.25
-pepo, wts, = time_evolve(pepo, ham, dt, nstep, alg, wts)
-env = converge_env(InfinitePartitionFunction(pepo), 16)
-result_2β = measure_mag(pepo, env)
+alg = SimpleUpdate(;
+    ψ0 = ψ, env0 = wts, H, dt, nstep,
+    trscheme = trscheme_ψ, gate_bothsides = true
+)
+ψ, wts, = time_evolve(alg)
+env = converge_env(InfinitePartitionFunction(ψ), 16)
+result_2β = measure_mag(ψ, env)
 @info "Magnetization at T = $(1 / (2β))" result_2β
 @test isapprox(abs.(result_2β), bm_2β, rtol = 1.0e-4)
 
 # Purification approach: results at 2β, or T = 1.25
-alg = SimpleUpdate(; trscheme = trscheme_pepo, gate_bothsides = false)
-pepo, = time_evolve(pepo0, ham, dt, 2 * nstep, alg, wts0)
-env = converge_env(InfinitePEPS(pepo), 8)
-result_2β′ = measure_mag(pepo, env; purified = true)
+alg = SimpleUpdate(;
+    ψ0, env0 = wts0, H, dt, nstep = 2 * nstep,
+    trscheme = trscheme_ψ, gate_bothsides = false
+)
+ψ, = time_evolve(alg)
+env = converge_env(InfinitePEPS(ψ), 8)
+result_2β′ = measure_mag(ψ, env; purified = true)
 @info "Magnetization at T = $(1 / (2β)) (purification approach)" result_2β′
 @test isapprox(abs.(result_2β′), bm_2β, rtol = 1.0e-2)
diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl
index 90220aa8..65a8b893 100644
--- a/test/timeevol/cluster_projectors.jl
+++ b/test/timeevol/cluster_projectors.jl
@@ -108,7 +108,7 @@ end
     trscheme_env = truncerr(1.0e-12) & truncdim(16)
     peps = InfinitePEPS(rand, Float64, Pspace, Vspace; unitcell = (Nr, Nc))
     wts = SUWeight(peps)
-    ham = real(
+    H = real(
         hubbard_model(
             ComplexF64, Trivial, U1Irrep, InfiniteSquare(Nr, Nc); t = 1.0, U = 8.0, mu = 0.0
         ),
@@ -119,25 +119,33 @@ end
     nstep = 5000
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
         trscheme = truncerr(1.0e-10) & truncdim(n == 1 ? 4 : 2)
-        alg = SimpleUpdate(; tol, trscheme, bipartite = true, check_interval = 1000)
-        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts)
+        alg = SimpleUpdate(;
+            ψ0 = peps, env0 = wts, H, dt, nstep, tol, trscheme,
+            bipartite = true, check_interval = 1000
+        )
+        peps, wts, = time_evolve(alg)
     end
     normalize!.(peps.A, Inf)
-    env = CTMRGEnv(rand, Float64, peps, Espace)
-    env, = leading_boundary(env, peps; tol = ctmrg_tol, trscheme = trscheme_env)
-    e_site = cost_function(peps, env, ham) / (Nr * Nc)
+    env, = leading_boundary(
+        CTMRGEnv(wts, peps), peps;
+        tol = ctmrg_tol, trscheme = trscheme_env
+    )
+    e_site = cost_function(peps, env, H) / (Nr * Nc)
     @info "2-site simple update energy = $e_site"
     # continue with 3-site simple update; energy should not change much
     dts = [1.0e-2, 5.0e-3]
     tols = [1.0e-8, 1.0e-8]
     trscheme = truncerr(1.0e-10) & truncdim(2)
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
-        alg = SimpleUpdate(; tol, trscheme, check_interval = 1000, force_3site = true)
-        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts)
+        alg = SimpleUpdate(;
+            ψ0 = peps, env0 = wts, H, dt, nstep, tol, trscheme,
+            check_interval = 1000, force_3site = true
+        )
+        peps, wts, = time_evolve(alg)
     end
     normalize!.(peps.A, Inf)
     env, = leading_boundary(env, peps; tol = ctmrg_tol, trscheme = trscheme_env)
-    e_site2 = cost_function(peps, env, ham) / (Nr * Nc)
+    e_site2 = cost_function(peps, env, H) / (Nr * Nc)
     @info "3-site simple update energy = $e_site2"
     @test e_site ≈ e_site2 atol = 5.0e-4
 end
diff --git a/test/timeevol/sitedep_truncation.jl b/test/timeevol/sitedep_truncation.jl
index 3ee588d6..cb038f1a 100644
--- a/test/timeevol/sitedep_truncation.jl
+++ b/test/timeevol/sitedep_truncation.jl
@@ -19,7 +19,7 @@ end
 
 @testset "Simple update: bipartite 2-site" begin
     Nr, Nc = 2, 2
-    ham = real(heisenberg_XYZ(InfiniteSquare(Nr, Nc); Jx = 1.0, Jy = 1.0, Jz = 1.0))
+    H = real(heisenberg_XYZ(InfiniteSquare(Nr, Nc); Jx = 1.0, Jy = 1.0, Jz = 1.0))
     Random.seed!(100)
     peps0 = InfinitePEPS(rand, Float64, ℂ^2, ℂ^10; unitcell = (Nr, Nc))
     env0 = SUWeight(peps0)
@@ -27,8 +27,8 @@ end
     # set trscheme to be compatible with bipartite structure
     bonddims = stack([[6 4; 4 6], [5 7; 7 5]]; dims = 1)
     trscheme = SiteDependentTruncation(collect(truncdim(d) for d in bonddims))
-    alg = SimpleUpdate(; trscheme, bipartite = true)
-    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
+    alg = SimpleUpdate(; ψ0 = peps0, env0, H, dt = 1.0e-2, nstep = 4, trscheme, bipartite = true)
+    peps, env, = time_evolve(alg)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
     # check bipartite structure is preserved
@@ -51,15 +51,16 @@ end
     # Site dependent truncation
     bonddims = rand(2:8, 2, Nr, Nc)
     trscheme = SiteDependentTruncation(collect(truncdim(d) for d in bonddims))
-    alg = SimpleUpdate(; trscheme)
     # 2-site SU
-    ham = real(heisenberg_XYZ(InfiniteSquare(Nr, Nc); Jx = 1.0, Jy = 1.0, Jz = 1.0))
-    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
+    H = real(heisenberg_XYZ(InfiniteSquare(Nr, Nc); Jx = 1.0, Jy = 1.0, Jz = 1.0))
+    alg = SimpleUpdate(; ψ0 = peps0, env0, H, dt = 1.0e-2, nstep = 4, trscheme)
+    peps, env, = time_evolve(alg)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
     # 3-site SU
-    ham = real(j1_j2_model(InfiniteSquare(Nr, Nc); J1 = 1.0, J2 = 0.5, sublattice = false))
-    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
+    H = real(j1_j2_model(InfiniteSquare(Nr, Nc); J1 = 1.0, J2 = 0.5, sublattice = false))
+    alg = SimpleUpdate(; ψ0 = peps0, env0, H, dt = 1.0e-2, nstep = 4, trscheme)
+    peps, env, = time_evolve(alg)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
 end

From b5468eed58a76e93a6fe82b70f364124208506dc Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 28 Oct 2025 22:59:17 +0800
Subject: [PATCH 06/36] Fix test

---
 test/bondenv/benv_fu.jl | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/test/bondenv/benv_fu.jl b/test/bondenv/benv_fu.jl
index 52bfca5c..d3e61530 100644
--- a/test/bondenv/benv_fu.jl
+++ b/test/bondenv/benv_fu.jl
@@ -13,8 +13,11 @@ function get_hubbard_state(t::Float64 = 1.0, U::Float64 = 8.0)
     Vphy = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 2, (1, 1 // 2) => 1, (1, -1 // 2) => 1)
     peps = InfinitePEPS(rand, ComplexF64, Vphy, Vphy; unitcell = (Nr, Nc))
     wts = SUWeight(peps)
-    alg = SimpleUpdate(; tol = 1.0e-8, trscheme = truncerr(1.0e-10) & truncdim(4), check_interval = 2000)
-    peps, = time_evolve(peps, H, 1.0e-2, 10000, alg, wts)
+    alg = SimpleUpdate(;
+        ψ0 = peps, env0 = wts, H, dt = 1.0e-2, nstep = 10000, tol = 1.0e-8,
+        trscheme = truncerr(1.0e-10) & truncdim(4), check_interval = 2000
+    )
+    peps, = time_evolve(alg)
     normalize!.(peps.A, Inf)
     return peps
 end

From f4b98a7ab145224009f66b5291a05c911c0ddcea Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Wed, 29 Oct 2025 16:59:54 +0800
Subject: [PATCH 07/36] Minor fixes [skip ci]

---
 src/algorithms/time_evolution/simpleupdate3site.jl | 2 +-
 src/algorithms/time_evolution/time_evolve.jl       | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl
index abd37e4b..cba5d51f 100644
--- a/src/algorithms/time_evolution/simpleupdate3site.jl
+++ b/src/algorithms/time_evolution/simpleupdate3site.jl
@@ -95,7 +95,7 @@ Here `-*-` is the twist on the physical axis.
 # Truncation of a bond on OBC-MPS
 
 Suppose we want to truncate the bond between 
-the n-th and the (n+1)-th sites such that the truncated ψ
+the n-th and the (n+1)-th sites such that the truncated state
 ```
     |ψ̃⟩  =  M[1]---...---M̃[n]---M̃[n+1]---...---M[N]
             ↓            ↓      ↓              ↓
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index 7809e7fd..6ef53339 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -7,6 +7,6 @@ function MPSKit.time_evolve(alg::Alg) where {Alg <: TimeEvolution}
         result = state
     end
     time_end = time()
-    @info @sprintf("Simple update finished. Total time elasped: %.2f s", time_end - time_start)
+    @info @sprintf("Time evolution finished. Total time elasped: %.2f s", time_end - time_start)
     return result
 end

From 8f7dff4594d25dd15564511b56db84a26f9ec6e9 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 4 Nov 2025 10:02:35 +0800
Subject: [PATCH 08/36] Fix merge regressions

---
 examples/j1j2_su/main.jl                      |   2 +-
 .../time_evolution/simpleupdate3site.jl       |   2 +-
 test/timeevol/cluster_projectors.jl           | 154 +++++++++---------
 test/timeevol/sitedep_truncation.jl           |   4 +-
 4 files changed, 81 insertions(+), 81 deletions(-)

diff --git a/examples/j1j2_su/main.jl b/examples/j1j2_su/main.jl
index d2e34814..d50e256e 100644
--- a/examples/j1j2_su/main.jl
+++ b/examples/j1j2_su/main.jl
@@ -119,7 +119,7 @@ using MPSKit: randomize!
 noise_peps = InfinitePEPS(randomize!.(deepcopy(peps.A)))
 peps₀ = peps + 1.0e-1noise_peps
 peps_opt, env_opt, E_opt, = fixedpoint(
-    H, peps₀, env; optimizer_alg = (; tol = 1.0e-4, nstep = 80)
+    H, peps₀, env; optimizer_alg = (; tol = 1.0e-4, maxiter = 80)
 );
 
 md"""
diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl
index 6e3ede02..67e0f38c 100644
--- a/src/algorithms/time_evolution/simpleupdate3site.jl
+++ b/src/algorithms/time_evolution/simpleupdate3site.jl
@@ -463,7 +463,7 @@ function _su3site_se!(
         ψ::InfiniteState, gs::Vector{T}, env::SUWeight,
         row::Int, col::Int, truncs::Vector{E};
         gate_bothsides::Bool = true
-    ) where {T <: AbstractTensorMap, E <: TruncationScheme}
+    ) where {T <: AbstractTensorMap, E <: TruncationStrategy}
     Nr, Nc = size(ψ)
     @assert 1 <= row <= Nr && 1 <= col <= Nc
     rm1, cp1 = _prev(row, Nr), _next(col, Nc)
diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl
index 1f2d0f6d..f710c4c2 100644
--- a/test/timeevol/cluster_projectors.jl
+++ b/test/timeevol/cluster_projectors.jl
@@ -21,82 +21,82 @@ Vspaces = [
     ),
 ]
 
-@testset "Cluster bond truncation with projectors" begin
-    Random.seed!(0)
-    N, n = 5, 2
-    for (Vphy, Vns, V) in Vspaces
-        Vvirs = fill(Vns, N + 1)
-        Vvirs[n + 1] = V
-        Ms1 = map(1:N) do i
-            Vw, Ve = Vvirs[i], Vvirs[i + 1]
-            return rand(Vw ⊗ Vphy ⊗ Vns' ⊗ Vns ← Ve)
-        end
-        normalize!.(Ms1, Inf)
-        revs = [isdual(space(M, 1)) for M in Ms1[2:end]]
-        # no truncation
-        Ms2 = deepcopy(Ms1)
-        wts2, ϵs, = _cluster_truncate!(Ms2, fill(FixedSpaceTruncation(), N - 1), revs)
-        @test all((ϵ == 0) for ϵ in ϵs)
-        normalize!.(Ms2, Inf)
-        @test fidelity_cluster(Ms1, Ms2) ≈ 1.0
-        lorths, rorths = verify_cluster_orth(Ms2, wts2)
-        @test all(lorths) && all(rorths)
-        # truncation on one bond
-        Ms3 = deepcopy(Ms1)
-        tspace = isdual(Vns) ? flip(Vns) : Vns
-        wts3, ϵs, = _cluster_truncate!(Ms3, fill(truncspace(tspace), N - 1), revs)
-        @test all((i == n) || (ϵ == 0) for (i, ϵ) in enumerate(ϵs))
-        normalize!.(Ms3, Inf)
-        ϵ = ϵs[n]
-        wt2, wt3 = wts2[n], wts3[n]
-        fid3, fid3_ = fidelity_cluster(Ms1, Ms3), fidelity_cluster(Ms2, Ms3)
-        @info "Fidelity of truncated cluster = $(fid3)"
-        @test fid3 ≈ fid3_
-        @test fid3 ≈ (norm(wt3) / norm(wt2))^2
-        @test fid3 ≈ 1.0 - (ϵ / norm(wt2))^2
-    end
-end
+# @testset "Cluster bond truncation with projectors" begin
+#     Random.seed!(0)
+#     N, n = 5, 2
+#     for (Vphy, Vns, V) in Vspaces
+#         Vvirs = fill(Vns, N + 1)
+#         Vvirs[n + 1] = V
+#         Ms1 = map(1:N) do i
+#             Vw, Ve = Vvirs[i], Vvirs[i + 1]
+#             return rand(Vw ⊗ Vphy ⊗ Vns' ⊗ Vns ← Ve)
+#         end
+#         normalize!.(Ms1, Inf)
+#         revs = [isdual(space(M, 1)) for M in Ms1[2:end]]
+#         # no truncation
+#         Ms2 = deepcopy(Ms1)
+#         wts2, ϵs, = _cluster_truncate!(Ms2, fill(FixedSpaceTruncation(), N - 1), revs)
+#         @test all((ϵ == 0) for ϵ in ϵs)
+#         normalize!.(Ms2, Inf)
+#         @test fidelity_cluster(Ms1, Ms2) ≈ 1.0
+#         lorths, rorths = verify_cluster_orth(Ms2, wts2)
+#         @test all(lorths) && all(rorths)
+#         # truncation on one bond
+#         Ms3 = deepcopy(Ms1)
+#         tspace = isdual(Vns) ? flip(Vns) : Vns
+#         wts3, ϵs, = _cluster_truncate!(Ms3, fill(truncspace(tspace), N - 1), revs)
+#         @test all((i == n) || (ϵ == 0) for (i, ϵ) in enumerate(ϵs))
+#         normalize!.(Ms3, Inf)
+#         ϵ = ϵs[n]
+#         wt2, wt3 = wts2[n], wts3[n]
+#         fid3, fid3_ = fidelity_cluster(Ms1, Ms3), fidelity_cluster(Ms2, Ms3)
+#         @info "Fidelity of truncated cluster = $(fid3)"
+#         @test fid3 ≈ fid3_
+#         @test fid3 ≈ (norm(wt3) / norm(wt2))^2
+#         @test fid3 ≈ 1.0 - (ϵ / norm(wt2))^2
+#     end
+# end
 
-@testset "Identity gate on 3-site cluster" begin
-    N, n = 3, 1
-    for (Vphy, Vns, V) in Vspaces
-        Vvirs = fill(Vns, N + 1)
-        Vvirs[n + 1] = V
-        Ms1 = map(1:N) do i
-            Vw, Ve = Vvirs[i], Vvirs[i + 1]
-            return normalize(rand(Vw ⊗ Vphy ⊗ Vns' ⊗ Vns ← Ve), Inf)
-        end
-        unit = id(Vphy)
-        gate = reduce(⊗, fill(unit, 3))
-        gs = PEPSKit.gate_to_mpo3(gate)
-        @test mpo_to_gate3(gs) ≈ gate
-        Ms2 = deepcopy(Ms1)
-        PEPSKit._apply_gatempo!(Ms2, gs)
-        fid = fidelity_cluster(Ms1, Ms2)
-        @test fid ≈ 1.0
-    end
-    for (Vphy, Vns, V) in Vspaces
-        Vvirs = fill(Vns, N + 1)
-        Vvirs[n + 1] = V
-        Ms1 = map(1:N) do i
-            Vw, Ve = Vvirs[i], Vvirs[i + 1]
-            return normalize(rand(Vw ⊗ Vphy ⊗ Vphy' ⊗ Vns' ⊗ Vns ← Ve), Inf)
-        end
-        unit = id(Vphy)
-        gate = reduce(⊗, fill(unit, 3))
-        gs = PEPSKit.gate_to_mpo3(gate)
-        @test mpo_to_gate3(gs) ≈ gate
-        for gate_ax in 1:2
-            Ms2 = deepcopy(Ms1)
-            PEPSKit._apply_gatempo!(Ms2, gs)
-            fid = fidelity_cluster(
-                [first(PEPSKit._fuse_physicalspaces(M)) for M in Ms1],
-                [first(PEPSKit._fuse_physicalspaces(M)) for M in Ms2]
-            )
-            @test fid ≈ 1.0
-        end
-    end
-end
+# @testset "Identity gate on 3-site cluster" begin
+#     N, n = 3, 1
+#     for (Vphy, Vns, V) in Vspaces
+#         Vvirs = fill(Vns, N + 1)
+#         Vvirs[n + 1] = V
+#         Ms1 = map(1:N) do i
+#             Vw, Ve = Vvirs[i], Vvirs[i + 1]
+#             return normalize(rand(Vw ⊗ Vphy ⊗ Vns' ⊗ Vns ← Ve), Inf)
+#         end
+#         unit = id(Vphy)
+#         gate = reduce(⊗, fill(unit, 3))
+#         gs = PEPSKit.gate_to_mpo3(gate)
+#         @test mpo_to_gate3(gs) ≈ gate
+#         Ms2 = deepcopy(Ms1)
+#         PEPSKit._apply_gatempo!(Ms2, gs)
+#         fid = fidelity_cluster(Ms1, Ms2)
+#         @test fid ≈ 1.0
+#     end
+#     for (Vphy, Vns, V) in Vspaces
+#         Vvirs = fill(Vns, N + 1)
+#         Vvirs[n + 1] = V
+#         Ms1 = map(1:N) do i
+#             Vw, Ve = Vvirs[i], Vvirs[i + 1]
+#             return normalize(rand(Vw ⊗ Vphy ⊗ Vphy' ⊗ Vns' ⊗ Vns ← Ve), Inf)
+#         end
+#         unit = id(Vphy)
+#         gate = reduce(⊗, fill(unit, 3))
+#         gs = PEPSKit.gate_to_mpo3(gate)
+#         @test mpo_to_gate3(gs) ≈ gate
+#         for gate_ax in 1:2
+#             Ms2 = deepcopy(Ms1)
+#             PEPSKit._apply_gatempo!(Ms2, gs)
+#             fid = fidelity_cluster(
+#                 [first(PEPSKit._fuse_physicalspaces(M)) for M in Ms1],
+#                 [first(PEPSKit._fuse_physicalspaces(M)) for M in Ms2]
+#             )
+#             @test fid ≈ 1.0
+#         end
+#     end
+# end
 
 @testset "Hubbard model with 2-site and 3-site SU" begin
     Nr, Nc = 2, 2
@@ -129,7 +129,7 @@ end
     normalize!.(peps.A, Inf)
     env, = leading_boundary(
         CTMRGEnv(wts, peps), peps;
-        tol = ctmrg_tol, trscheme = trscheme_env
+        tol = ctmrg_tol, trunc = trunc_env
     )
     e_site = cost_function(peps, env, H) / (Nr * Nc)
     @info "2-site simple update energy = $e_site"
@@ -145,7 +145,7 @@ end
         peps, wts, = time_evolve(alg)
     end
     normalize!.(peps.A, Inf)
-    env, = leading_boundary(env, peps; tol = ctmrg_tol, trscheme = trscheme_env)
+    env, = leading_boundary(env, peps; tol = ctmrg_tol, trunc = trunc_env)
     e_site2 = cost_function(peps, env, H) / (Nr * Nc)
     @info "3-site simple update energy = $e_site2"
     @test e_site ≈ e_site2 atol = 5.0e-4
diff --git a/test/timeevol/sitedep_truncation.jl b/test/timeevol/sitedep_truncation.jl
index fe9f97d1..3f7a48ad 100644
--- a/test/timeevol/sitedep_truncation.jl
+++ b/test/timeevol/sitedep_truncation.jl
@@ -54,13 +54,13 @@ end
     trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims))
     # 2-site SU
     H = real(heisenberg_XYZ(InfiniteSquare(Nr, Nc); Jx = 1.0, Jy = 1.0, Jz = 1.0))
-    alg = SimpleUpdate(; ψ0 = peps0, env0, H, dt = 1.0e-2, nstep = 4, trscheme)
+    alg = SimpleUpdate(; ψ0 = peps0, env0, H, dt = 1.0e-2, nstep = 4, trunc)
     peps, env, = time_evolve(alg)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
     # 3-site SU
     H = real(j1_j2_model(InfiniteSquare(Nr, Nc); J1 = 1.0, J2 = 0.5, sublattice = false))
-    alg = SimpleUpdate(; ψ0 = peps0, env0, H, dt = 1.0e-2, nstep = 4, trscheme)
+    alg = SimpleUpdate(; ψ0 = peps0, env0, H, dt = 1.0e-2, nstep = 4, trunc)
     peps, env, = time_evolve(alg)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims

From 542b58efeef54fc871486861633f9c57dd94a707 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 4 Nov 2025 10:29:23 +0800
Subject: [PATCH 09/36] Restore accidentally commented tests

---
 test/timeevol/cluster_projectors.jl | 150 ++++++++++++++--------------
 1 file changed, 75 insertions(+), 75 deletions(-)

diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl
index f710c4c2..9a091cca 100644
--- a/test/timeevol/cluster_projectors.jl
+++ b/test/timeevol/cluster_projectors.jl
@@ -21,82 +21,82 @@ Vspaces = [
     ),
 ]
 
-# @testset "Cluster bond truncation with projectors" begin
-#     Random.seed!(0)
-#     N, n = 5, 2
-#     for (Vphy, Vns, V) in Vspaces
-#         Vvirs = fill(Vns, N + 1)
-#         Vvirs[n + 1] = V
-#         Ms1 = map(1:N) do i
-#             Vw, Ve = Vvirs[i], Vvirs[i + 1]
-#             return rand(Vw ⊗ Vphy ⊗ Vns' ⊗ Vns ← Ve)
-#         end
-#         normalize!.(Ms1, Inf)
-#         revs = [isdual(space(M, 1)) for M in Ms1[2:end]]
-#         # no truncation
-#         Ms2 = deepcopy(Ms1)
-#         wts2, ϵs, = _cluster_truncate!(Ms2, fill(FixedSpaceTruncation(), N - 1), revs)
-#         @test all((ϵ == 0) for ϵ in ϵs)
-#         normalize!.(Ms2, Inf)
-#         @test fidelity_cluster(Ms1, Ms2) ≈ 1.0
-#         lorths, rorths = verify_cluster_orth(Ms2, wts2)
-#         @test all(lorths) && all(rorths)
-#         # truncation on one bond
-#         Ms3 = deepcopy(Ms1)
-#         tspace = isdual(Vns) ? flip(Vns) : Vns
-#         wts3, ϵs, = _cluster_truncate!(Ms3, fill(truncspace(tspace), N - 1), revs)
-#         @test all((i == n) || (ϵ == 0) for (i, ϵ) in enumerate(ϵs))
-#         normalize!.(Ms3, Inf)
-#         ϵ = ϵs[n]
-#         wt2, wt3 = wts2[n], wts3[n]
-#         fid3, fid3_ = fidelity_cluster(Ms1, Ms3), fidelity_cluster(Ms2, Ms3)
-#         @info "Fidelity of truncated cluster = $(fid3)"
-#         @test fid3 ≈ fid3_
-#         @test fid3 ≈ (norm(wt3) / norm(wt2))^2
-#         @test fid3 ≈ 1.0 - (ϵ / norm(wt2))^2
-#     end
-# end
+@testset "Cluster bond truncation with projectors" begin
+    Random.seed!(0)
+    N, n = 5, 2
+    for (Vphy, Vns, V) in Vspaces
+        Vvirs = fill(Vns, N + 1)
+        Vvirs[n + 1] = V
+        Ms1 = map(1:N) do i
+            Vw, Ve = Vvirs[i], Vvirs[i + 1]
+            return rand(Vw ⊗ Vphy ⊗ Vns' ⊗ Vns ← Ve)
+        end
+        normalize!.(Ms1, Inf)
+        revs = [isdual(space(M, 1)) for M in Ms1[2:end]]
+        # no truncation
+        Ms2 = deepcopy(Ms1)
+        wts2, ϵs, = _cluster_truncate!(Ms2, fill(FixedSpaceTruncation(), N - 1), revs)
+        @test all((ϵ == 0) for ϵ in ϵs)
+        normalize!.(Ms2, Inf)
+        @test fidelity_cluster(Ms1, Ms2) ≈ 1.0
+        lorths, rorths = verify_cluster_orth(Ms2, wts2)
+        @test all(lorths) && all(rorths)
+        # truncation on one bond
+        Ms3 = deepcopy(Ms1)
+        tspace = isdual(Vns) ? flip(Vns) : Vns
+        wts3, ϵs, = _cluster_truncate!(Ms3, fill(truncspace(tspace), N - 1), revs)
+        @test all((i == n) || (ϵ == 0) for (i, ϵ) in enumerate(ϵs))
+        normalize!.(Ms3, Inf)
+        ϵ = ϵs[n]
+        wt2, wt3 = wts2[n], wts3[n]
+        fid3, fid3_ = fidelity_cluster(Ms1, Ms3), fidelity_cluster(Ms2, Ms3)
+        @info "Fidelity of truncated cluster = $(fid3)"
+        @test fid3 ≈ fid3_
+        @test fid3 ≈ (norm(wt3) / norm(wt2))^2
+        @test fid3 ≈ 1.0 - (ϵ / norm(wt2))^2
+    end
+end
 
-# @testset "Identity gate on 3-site cluster" begin
-#     N, n = 3, 1
-#     for (Vphy, Vns, V) in Vspaces
-#         Vvirs = fill(Vns, N + 1)
-#         Vvirs[n + 1] = V
-#         Ms1 = map(1:N) do i
-#             Vw, Ve = Vvirs[i], Vvirs[i + 1]
-#             return normalize(rand(Vw ⊗ Vphy ⊗ Vns' ⊗ Vns ← Ve), Inf)
-#         end
-#         unit = id(Vphy)
-#         gate = reduce(⊗, fill(unit, 3))
-#         gs = PEPSKit.gate_to_mpo3(gate)
-#         @test mpo_to_gate3(gs) ≈ gate
-#         Ms2 = deepcopy(Ms1)
-#         PEPSKit._apply_gatempo!(Ms2, gs)
-#         fid = fidelity_cluster(Ms1, Ms2)
-#         @test fid ≈ 1.0
-#     end
-#     for (Vphy, Vns, V) in Vspaces
-#         Vvirs = fill(Vns, N + 1)
-#         Vvirs[n + 1] = V
-#         Ms1 = map(1:N) do i
-#             Vw, Ve = Vvirs[i], Vvirs[i + 1]
-#             return normalize(rand(Vw ⊗ Vphy ⊗ Vphy' ⊗ Vns' ⊗ Vns ← Ve), Inf)
-#         end
-#         unit = id(Vphy)
-#         gate = reduce(⊗, fill(unit, 3))
-#         gs = PEPSKit.gate_to_mpo3(gate)
-#         @test mpo_to_gate3(gs) ≈ gate
-#         for gate_ax in 1:2
-#             Ms2 = deepcopy(Ms1)
-#             PEPSKit._apply_gatempo!(Ms2, gs)
-#             fid = fidelity_cluster(
-#                 [first(PEPSKit._fuse_physicalspaces(M)) for M in Ms1],
-#                 [first(PEPSKit._fuse_physicalspaces(M)) for M in Ms2]
-#             )
-#             @test fid ≈ 1.0
-#         end
-#     end
-# end
+@testset "Identity gate on 3-site cluster" begin
+    N, n = 3, 1
+    for (Vphy, Vns, V) in Vspaces
+        Vvirs = fill(Vns, N + 1)
+        Vvirs[n + 1] = V
+        Ms1 = map(1:N) do i
+            Vw, Ve = Vvirs[i], Vvirs[i + 1]
+            return normalize(rand(Vw ⊗ Vphy ⊗ Vns' ⊗ Vns ← Ve), Inf)
+        end
+        unit = id(Vphy)
+        gate = reduce(⊗, fill(unit, 3))
+        gs = PEPSKit.gate_to_mpo3(gate)
+        @test mpo_to_gate3(gs) ≈ gate
+        Ms2 = deepcopy(Ms1)
+        PEPSKit._apply_gatempo!(Ms2, gs)
+        fid = fidelity_cluster(Ms1, Ms2)
+        @test fid ≈ 1.0
+    end
+    for (Vphy, Vns, V) in Vspaces
+        Vvirs = fill(Vns, N + 1)
+        Vvirs[n + 1] = V
+        Ms1 = map(1:N) do i
+            Vw, Ve = Vvirs[i], Vvirs[i + 1]
+            return normalize(rand(Vw ⊗ Vphy ⊗ Vphy' ⊗ Vns' ⊗ Vns ← Ve), Inf)
+        end
+        unit = id(Vphy)
+        gate = reduce(⊗, fill(unit, 3))
+        gs = PEPSKit.gate_to_mpo3(gate)
+        @test mpo_to_gate3(gs) ≈ gate
+        for gate_ax in 1:2
+            Ms2 = deepcopy(Ms1)
+            PEPSKit._apply_gatempo!(Ms2, gs)
+            fid = fidelity_cluster(
+                [first(PEPSKit._fuse_physicalspaces(M)) for M in Ms1],
+                [first(PEPSKit._fuse_physicalspaces(M)) for M in Ms2]
+            )
+            @test fid ≈ 1.0
+        end
+    end
+end
 
 @testset "Hubbard model with 2-site and 3-site SU" begin
     Nr, Nc = 2, 2

From 6e65759696921eb0aa9b7fd5a1f0c24653ae1547 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 4 Nov 2025 15:13:00 +0800
Subject: [PATCH 10/36] Undo unnecessary renames

---
 src/algorithms/time_evolution/simpleupdate.jl | 38 +++++++-------
 .../time_evolution/simpleupdate3site.jl       | 40 +++++++--------
 test/examples/heisenberg.jl                   | 14 +++---
 test/examples/j1j2_finiteT.jl                 | 28 +++++------
 test/examples/tf_ising_finiteT.jl             | 50 +++++++++----------
 test/timeevol/cluster_projectors.jl           | 16 +++---
 test/timeevol/sitedep_truncation.jl           | 11 ++--
 7 files changed, 98 insertions(+), 99 deletions(-)

diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index 0461d40f..be56d442 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -45,14 +45,14 @@ Simple update of the x-bond between `[r,c]` and `[r,c+1]`.
 ```
 =#
 function _su_xbond!(
-        ψ::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
+        state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
         row::Int, col::Int, trunc::TruncationStrategy; gate_ax::Int = 1
     ) where {T <: Number, S <: ElementarySpace}
-    Nr, Nc = size(ψ)[1:2]
+    Nr, Nc = size(state)[1:2]
     @assert 1 <= row <= Nr && 1 <= col <= Nc
     cp1 = _next(col, Nc)
     # absorb environment weights
-    A, B = ψ.A[row, col], ψ.A[row, cp1]
+    A, B = state.A[row, col], state.A[row, cp1]
     A = absorb_weight(A, env, row, col, (NORTH, SOUTH, WEST); inv = false)
     B = absorb_weight(B, env, row, cp1, (NORTH, SOUTH, EAST); inv = false)
     normalize!(A, Inf)
@@ -68,7 +68,7 @@ function _su_xbond!(
     normalize!(B, Inf)
     normalize!(s, Inf)
     # update tensor dict and weight on current bond
-    ψ.A[row, col], ψ.A[row, cp1] = A, B
+    state.A[row, col], state.A[row, cp1] = A, B
     env.data[1, row, col] = s
     return ϵ
 end
@@ -84,14 +84,14 @@ Simple update of the y-bond between `[r,c]` and `[r-1,c]`.
 ```
 =#
 function _su_ybond!(
-        ψ::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
+        state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
         row::Int, col::Int, trunc::TruncationStrategy; gate_ax::Int = 1
     ) where {T <: Number, S <: ElementarySpace}
-    Nr, Nc = size(ψ)[1:2]
+    Nr, Nc = size(state)[1:2]
     @assert 1 <= row <= Nr && 1 <= col <= Nc
     rm1 = _prev(row, Nr)
     # absorb environment weights
-    A, B = ψ.A[row, col], ψ.A[rm1, col]
+    A, B = state.A[row, col], state.A[rm1, col]
     A = absorb_weight(A, env, row, col, (EAST, SOUTH, WEST); inv = false)
     B = absorb_weight(B, env, rm1, col, (NORTH, EAST, WEST); inv = false)
     normalize!(A, Inf)
@@ -107,48 +107,48 @@ function _su_ybond!(
     normalize!(A, Inf)
     normalize!(B, Inf)
     normalize!(s, Inf)
-    ψ.A[row, col], ψ.A[rm1, col] = A, B
+    state.A[row, col], state.A[rm1, col] = A, B
     env.data[2, row, col] = s
     return ϵ
 end
 
 function su_iter(
-        ψ::InfiniteState, gate::LocalOperator, alg::SimpleUpdate, env::SUWeight
+        state::InfiniteState, gate::LocalOperator, alg::SimpleUpdate, env::SUWeight
     )
-    @assert size(gate.lattice) == size(ψ)[1:2]
-    Nr, Nc = size(ψ)[1:2]
+    @assert size(gate.lattice) == size(state)[1:2]
+    Nr, Nc = size(state)[1:2]
     alg.bipartite && (@assert Nr == Nc == 2)
     (Nr >= 2 && Nc >= 2) || throw(
         ArgumentError("`state` unit cell size for simple update should be no smaller than (2, 2)."),
     )
-    ψ2, env2 = deepcopy(ψ), deepcopy(env)
+    state2, env2 = deepcopy(state), deepcopy(env)
     gate_axs = alg.gate_bothsides ? (1:2) : (1:1)
     for r in 1:Nr, c in 1:Nc
         term = get_gateterm(gate, (CartesianIndex(r, c), CartesianIndex(r, c + 1)))
         trunc = truncation_strategy(alg.trunc, 1, r, c)
         for gate_ax in gate_axs
-            _su_xbond!(ψ2, term, env2, r, c, trunc; gate_ax)
+            _su_xbond!(state2, term, env2, r, c, trunc; gate_ax)
         end
         if alg.bipartite
             rp1, cp1 = _next(r, Nr), _next(c, Nc)
-            ψ2.A[rp1, cp1] = deepcopy(ψ2.A[r, c])
-            ψ2.A[rp1, c] = deepcopy(ψ2.A[r, cp1])
+            state2.A[rp1, cp1] = deepcopy(state2.A[r, c])
+            state2.A[rp1, c] = deepcopy(state2.A[r, cp1])
             env2.data[1, rp1, cp1] = deepcopy(env2.data[1, r, c])
         end
         term = get_gateterm(gate, (CartesianIndex(r, c), CartesianIndex(r - 1, c)))
         trunc = truncation_strategy(alg.trunc, 2, r, c)
         for gate_ax in gate_axs
-            _su_ybond!(ψ2, term, env2, r, c, trunc; gate_ax)
+            _su_ybond!(state2, term, env2, r, c, trunc; gate_ax)
         end
         if alg.bipartite
             rm1, cm1 = _prev(r, Nr), _prev(c, Nc)
-            ψ2.A[rm1, cm1] = deepcopy(ψ2.A[r, c])
-            ψ2.A[r, cm1] = deepcopy(ψ2.A[rm1, c])
+            state2.A[rm1, cm1] = deepcopy(state2.A[r, c])
+            state2.A[r, cm1] = deepcopy(state2.A[rm1, c])
             env2.data[2, rm1, cm1] = deepcopy(env2.data[2, r, c])
         end
     end
     wtdiff = compare_weights(env2, env)
-    return ψ2, env2, (; wtdiff)
+    return state2, env2, (; wtdiff)
 end
 
 # Initial iteration
diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl
index 67e0f38c..7dcbdd08 100644
--- a/src/algorithms/time_evolution/simpleupdate3site.jl
+++ b/src/algorithms/time_evolution/simpleupdate3site.jl
@@ -447,28 +447,28 @@ Obtain the 3-site cluster in the "southeast corner" of a square plaquette.
         c      c+1
 ```
 =#
-function get_3site_se(ψ::InfiniteState, env::SUWeight, row::Int, col::Int)
-    Nr, Nc = size(ψ)
+function get_3site_se(state::InfiniteState, env::SUWeight, row::Int, col::Int)
+    Nr, Nc = size(state)
     rm1, cp1 = _prev(row, Nr), _next(col, Nc)
     coords_se = [(row, col), (row, cp1), (rm1, cp1)]
-    perms_se = isa(ψ, InfinitePEPS) ? perms_se_peps : perms_se_pepo
+    perms_se = isa(state, InfinitePEPS) ? perms_se_peps : perms_se_pepo
     Ms = map(zip(coords_se, perms_se, openaxs_se)) do (coord, perm, openaxs)
-        M = absorb_weight(ψ.A[CartesianIndex(coord)], env, coord[1], coord[2], openaxs)
+        M = absorb_weight(state.A[CartesianIndex(coord)], env, coord[1], coord[2], openaxs)
         return permute(M, perm)
     end
     return Ms
 end
 
 function _su3site_se!(
-        ψ::InfiniteState, gs::Vector{T}, env::SUWeight,
+        state::InfiniteState, gs::Vector{T}, env::SUWeight,
         row::Int, col::Int, truncs::Vector{E};
         gate_bothsides::Bool = true
     ) where {T <: AbstractTensorMap, E <: TruncationStrategy}
-    Nr, Nc = size(ψ)
+    Nr, Nc = size(state)
     @assert 1 <= row <= Nr && 1 <= col <= Nc
     rm1, cp1 = _prev(row, Nr), _next(col, Nc)
     # southwest 3-site cluster and arrow direction within it
-    Ms = get_3site_se(ψ, env, row, col)
+    Ms = get_3site_se(state, env, row, col)
     revs = [isdual(space(M, 1)) for M in Ms[2:end]]
     Vphys = [codomain(M, 2) for M in Ms]
     normalize!.(Ms, Inf)
@@ -481,55 +481,55 @@ function _su3site_se!(
     ϵs = nothing
     for gate_ax in gate_axs
         _apply_gatempo!(Ms, gs; gate_ax)
-        if isa(ψ, InfinitePEPO)
+        if isa(state, InfinitePEPO)
             Ms = [first(_fuse_physicalspaces(M)) for M in Ms]
         end
         wts, ϵs, = _cluster_truncate!(Ms, truncs, revs)
-        if isa(ψ, InfinitePEPO)
+        if isa(state, InfinitePEPO)
             Ms = [first(_unfuse_physicalspace(M, Vphy)) for (M, Vphy) in zip(Ms, Vphys)]
         end
         for (wt, wt_idx) in zip(wts, wt_idxs)
             env[CartesianIndex(wt_idx)] = normalize(wt, Inf)
         end
     end
-    invperms_se = isa(ψ, InfinitePEPS) ? invperms_se_peps : invperms_se_pepo
+    invperms_se = isa(state, InfinitePEPS) ? invperms_se_peps : invperms_se_pepo
     for (M, coord, invperm, openaxs, Vphy) in zip(Ms, coords, invperms_se, openaxs_se, Vphys)
         # restore original axes order
         M = permute(M, invperm)
         # remove weights on open axes of the cluster
         M = absorb_weight(M, env, coord[1], coord[2], openaxs; inv = true)
-        ψ.A[CartesianIndex(coord)] = normalize(M, Inf)
+        state.A[CartesianIndex(coord)] = normalize(M, Inf)
     end
     return ϵs
 end
 
 function su_iter(
-        ψ::InfiniteState, gatempos::Vector{G}, alg::SimpleUpdate, env::SUWeight
+        state::InfiniteState, gatempos::Vector{G}, alg::SimpleUpdate, env::SUWeight
     ) where {G <: AbstractMatrix}
-    if ψ isa InfinitePEPO
-        @assert size(ψ, 3) == 1
+    if state isa InfinitePEPO
+        @assert size(state, 3) == 1
     end
-    Nr, Nc = size(ψ)[1:2]
+    Nr, Nc = size(state)[1:2]
     (Nr >= 2 && Nc >= 2) || throw(
         ArgumentError(
             "iPEPS unit cell size for simple update should be no smaller than (2, 2)."
         ),
     )
-    ψ2, env2 = deepcopy(ψ), deepcopy(env)
+    state2, env2 = deepcopy(state), deepcopy(env)
     trunc = alg.trunc
     for i in 1:4
-        Nr, Nc = size(ψ2)[1:2]
+        Nr, Nc = size(state2)[1:2]
         for r in 1:Nr, c in 1:Nc
             gs = gatempos[i][r, c]
             truncs = [
                 truncation_strategy(trunc, 1, r, c)
                 truncation_strategy(trunc, 2, r, _next(c, Nc))
             ]
-            _su3site_se!(ψ2, gs, env2, r, c, truncs; alg.gate_bothsides)
+            _su3site_se!(state2, gs, env2, r, c, truncs; alg.gate_bothsides)
         end
-        ψ2, env2 = rotl90(ψ2), rotl90(env2)
+        state2, env2 = rotl90(state2), rotl90(env2)
         trunc = rotl90(trunc)
     end
     wtdiff = compare_weights(env2, env)
-    return ψ2, env2, (; wtdiff)
+    return state2, env2, (; wtdiff)
 end
diff --git a/test/examples/heisenberg.jl b/test/examples/heisenberg.jl
index 3762c0e8..d4224af9 100644
--- a/test/examples/heisenberg.jl
+++ b/test/examples/heisenberg.jl
@@ -59,10 +59,10 @@ end
     normalize!.(peps.A, Inf)
 
     # Heisenberg model Hamiltonian
-    H = heisenberg_XYZ(InfiniteSquare(N1, N2); Jx = 1.0, Jy = 1.0, Jz = 1.0)
+    ham = heisenberg_XYZ(InfiniteSquare(N1, N2); Jx = 1.0, Jy = 1.0, Jz = 1.0)
     # assert imaginary part is zero
-    @assert length(imag(H).terms) == 0
-    H = real(H)
+    @assert length(imag(ham).terms) == 0
+    ham = real(ham)
 
     # simple update
     dts = [1.0e-2, 1.0e-3, 1.0e-3, 1.0e-4]
@@ -71,21 +71,23 @@ end
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
         Dbond2 = (n == 2) ? Dbond + 2 : Dbond
         trunc = truncerror(; atol = 1.0e-10) & truncrank(Dbond2)
-        alg = SimpleUpdate(; ψ0 = peps, env0 = wts, H, dt, nstep, trunc, bipartite = false)
+        alg = SimpleUpdate(; ψ0 = peps, env0 = wts, H = ham, dt, nstep, trunc, bipartite = false)
         peps, wts, = time_evolve(alg)
     end
 
     # measure physical quantities with CTMRG
     normalize!.(peps.A, Inf)
     env, = leading_boundary(CTMRGEnv(rand, Float64, peps, Espace), peps; tol = ctmrg_tol, maxiter = ctmrg_maxiter)
-    e_site = cost_function(peps, env, H) / (N1 * N2)
+    e_site = cost_function(peps, env, ham) / (N1 * N2)
     @info "Simple update energy = $e_site"
     # benchmark data from Phys. Rev. B 94, 035133 (2016)
     @test isapprox(e_site, -0.6594; atol = 1.0e-3)
 
     # continue with auto differentiation
     peps_final, env_final, E_final, = fixedpoint(
-        H, peps, env;
+        ham,
+        peps,
+        env;
         optimizer_alg = (; tol = gradtol, maxiter = 25),
         boundary_alg = (; maxiter = ctmrg_maxiter),
         gradient_alg = (; alg = :linsolver, solver_alg = (; alg = :gmres)),
diff --git a/test/examples/j1j2_finiteT.jl b/test/examples/j1j2_finiteT.jl
index c12273e1..9a3828a8 100644
--- a/test/examples/j1j2_finiteT.jl
+++ b/test/examples/j1j2_finiteT.jl
@@ -20,12 +20,12 @@ function converge_env(state, χ::Int)
 end
 
 Nr, Nc = 2, 2
-H = j1_j2_model(
+ham = j1_j2_model(
     Float64, SU2Irrep, InfiniteSquare(Nr, Nc);
     J1 = 1.0, J2 = 0.5, sublattice = false
 )
-ψ0 = PEPSKit.infinite_temperature_density_matrix(H)
-wts0 = SUWeight(ψ0)
+pepo0 = PEPSKit.infinite_temperature_density_matrix(ham)
+wts0 = SUWeight(pepo0)
 # 7 = 1 (spin-0) + 2 x 3 (spin-1)
 trunc_pepo = truncrank(7) & truncerror(; atol = 1.0e-12)
 check_interval = 100
@@ -33,27 +33,27 @@ check_interval = 100
 # PEPO approach
 dt, nstep = 1.0e-3, 600
 alg = SimpleUpdate(;
-    ψ0, env0 = wts0, H, dt, nstep, trunc = trunc_pepo,
+    ψ0 = pepo0, env0 = wts0, H = ham, dt, nstep, trunc = trunc_pepo,
     gate_bothsides = true, check_interval
 )
-ψ, wts, = time_evolve(alg)
-env = converge_env(InfinitePartitionFunction(ψ), 16)
-energy = expectation_value(ψ, H, env) / (Nr * Nc)
+pepo, wts, = time_evolve(alg)
+env = converge_env(InfinitePartitionFunction(pepo), 16)
+energy = expectation_value(pepo, ham, env) / (Nr * Nc)
 @info "β = $(dt * nstep): tr(ρH) = $(energy)"
 @test energy ≈ bm[2] atol = 5.0e-3
 
 # PEPS (purified PEPO) approach
 alg = SimpleUpdate(;
-    ψ0, env0 = wts0, H, dt, nstep, trunc = trunc_pepo,
+    ψ0 = pepo0, env0 = wts0, H = ham, dt, nstep, trunc = trunc_pepo,
     gate_bothsides = false, check_interval
 )
-ψ, wts, = time_evolve(alg)
-env = converge_env(InfinitePartitionFunction(ψ), 16)
-energy = expectation_value(ψ, H, env) / (Nr * Nc)
+pepo, wts, = time_evolve(alg)
+env = converge_env(InfinitePartitionFunction(pepo), 16)
+energy = expectation_value(pepo, ham, env) / (Nr * Nc)
 @info "β = $(dt * nstep) / 2: tr(ρH) = $(energy)"
 @test energy ≈ bm[1] atol = 5.0e-3
 
-env = converge_env(InfinitePEPS(ψ), 16)
-energy = expectation_value(ψ, H, ψ, env) / (Nr * Nc)
-@info "β = $(dt * nstep): ⟨ρ|H|ρ⟩ = $(energy)"
+env = converge_env(InfinitePEPS(pepo), 16)
+energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc)
+@info "β = $(dt * nstep): ⟨ρ|ham|ρ⟩ = $(energy)"
 @test energy ≈ bm[2] atol = 5.0e-3
diff --git a/test/examples/tf_ising_finiteT.jl b/test/examples/tf_ising_finiteT.jl
index 6b5d3c50..f3217b90 100644
--- a/test/examples/tf_ising_finiteT.jl
+++ b/test/examples/tf_ising_finiteT.jl
@@ -23,34 +23,34 @@ function tfising_model(T::Type{<:Number}, lattice::InfiniteSquare; J = 1.0, g =
     return PEPSKit.nearest_neighbour_hamiltonian(spaces, gate)
 end
 
-function converge_env(ψ, χ::Int)
+function converge_env(state, χ::Int)
     trunc1 = truncrank(4) & truncerror(; atol = 1.0e-12)
-    env0 = CTMRGEnv(randn, Float64, ψ, ℂ^4)
-    env, = leading_boundary(env0, ψ; alg = :sequential, trunc = trunc1, tol = 1.0e-10)
+    env0 = CTMRGEnv(randn, Float64, state, ℂ^4)
+    env, = leading_boundary(env0, state; alg = :sequential, trunc = trunc1, tol = 1.0e-10)
     trunc2 = truncrank(χ) & truncerror(; atol = 1.0e-12)
-    env, = leading_boundary(env, ψ; alg = :sequential, trunc = trunc2, tol = 1.0e-10)
+    env, = leading_boundary(env, state; alg = :sequential, trunc = trunc2, tol = 1.0e-10)
     return env
 end
 
-function measure_mag(ψ::InfinitePEPO, env::CTMRGEnv; purified::Bool = false)
+function measure_mag(pepo::InfinitePEPO, env::CTMRGEnv; purified::Bool = false)
     r, c = 1, 1
-    lattice = physicalspace(ψ)
+    lattice = physicalspace(pepo)
     Mx = LocalOperator(lattice, ((r, c),) => σˣ(Float64, Trivial))
     Mz = LocalOperator(lattice, ((r, c),) => σᶻ(Float64, Trivial))
     if purified
-        magx = expectation_value(ψ, Mx, ψ, env)
-        magz = expectation_value(ψ, Mz, ψ, env)
+        magx = expectation_value(pepo, Mx, pepo, env)
+        magz = expectation_value(pepo, Mz, pepo, env)
     else
-        magx = expectation_value(ψ, Mx, env)
-        magz = expectation_value(ψ, Mz, env)
+        magx = expectation_value(pepo, Mx, env)
+        magz = expectation_value(pepo, Mz, env)
     end
     return [magx, magz]
 end
 
 Nr, Nc = 2, 2
-H = tfising_model(Float64, InfiniteSquare(Nr, Nc); J = 1.0, g = 2.0)
-ψ0 = PEPSKit.infinite_temperature_density_matrix(H)
-wts0 = SUWeight(ψ0)
+ham = tfising_model(Float64, InfiniteSquare(Nr, Nc); J = 1.0, g = 2.0)
+pepo0 = PEPSKit.infinite_temperature_density_matrix(ham)
+wts0 = SUWeight(pepo0)
 
 trunc_pepo = truncrank(8) & truncerror(; atol = 1.0e-12)
 
@@ -61,33 +61,33 @@ dt, nstep = 1.0e-3, 400
 
 # PEPO approach: results at β, or T = 2.5
 alg = SimpleUpdate(;
-    ψ0, env0 = wts0, H, dt, nstep,
+    ψ0 = pepo0, env0 = wts0, H = ham, dt, nstep,
     trunc = trunc_pepo, gate_bothsides = true
 )
-ψ, wts, = time_evolve(alg)
-env = converge_env(InfinitePartitionFunction(ψ), 16)
-result_β = measure_mag(ψ, env)
+pepo, wts, = time_evolve(alg)
+env = converge_env(InfinitePartitionFunction(pepo), 16)
+result_β = measure_mag(pepo, env)
 @info "Magnetization at T = $(1 / β)" result_β
 @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2)
 
 # continue to get results at 2β, or T = 1.25
 alg = SimpleUpdate(;
-    ψ0 = ψ, env0 = wts, H, dt, nstep,
+    ψ0 = pepo, env0 = wts, H = ham, dt, nstep,
     trunc = trunc_pepo, gate_bothsides = true
 )
-ψ, wts, = time_evolve(alg)
-env = converge_env(InfinitePartitionFunction(ψ), 16)
-result_2β = measure_mag(ψ, env)
+pepo, wts, = time_evolve(alg)
+env = converge_env(InfinitePartitionFunction(pepo), 16)
+result_2β = measure_mag(pepo, env)
 @info "Magnetization at T = $(1 / (2β))" result_2β
 @test isapprox(abs.(result_2β), bm_2β, rtol = 1.0e-4)
 
 # Purification approach: results at 2β, or T = 1.25
 alg = SimpleUpdate(;
-    ψ0, env0 = wts0, H, dt, nstep = 2 * nstep,
+    ψ0 = pepo0, env0 = wts0, H = ham, dt, nstep = 2 * nstep,
     trunc = trunc_pepo, gate_bothsides = false
 )
-ψ, = time_evolve(alg)
-env = converge_env(InfinitePEPS(ψ), 8)
-result_2β′ = measure_mag(ψ, env; purified = true)
+pepo, = time_evolve(alg)
+env = converge_env(InfinitePEPS(pepo), 8)
+result_2β′ = measure_mag(pepo, env; purified = true)
 @info "Magnetization at T = $(1 / (2β)) (purification approach)" result_2β′
 @test isapprox(abs.(result_2β′), bm_2β, rtol = 1.0e-2)
diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl
index 9a091cca..501e6334 100644
--- a/test/timeevol/cluster_projectors.jl
+++ b/test/timeevol/cluster_projectors.jl
@@ -109,7 +109,7 @@ end
     trunc_env = truncerror(; atol = 1.0e-12) & truncrank(16)
     peps = InfinitePEPS(rand, Float64, Pspace, Vspace; unitcell = (Nr, Nc))
     wts = SUWeight(peps)
-    H = real(
+    ham = real(
         hubbard_model(
             ComplexF64, Trivial, U1Irrep, InfiniteSquare(Nr, Nc); t = 1.0, U = 8.0, mu = 0.0
         ),
@@ -121,17 +121,15 @@ end
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
         trunc = truncerror(; atol = 1.0e-10) & truncrank(n == 1 ? 4 : 2)
         alg = SimpleUpdate(;
-            ψ0 = peps, env0 = wts, H, dt, nstep, tol, trunc,
+            ψ0 = peps, env0 = wts, H = ham, dt, nstep, tol, trunc,
             bipartite = true, check_interval = 1000
         )
         peps, wts, = time_evolve(alg)
     end
     normalize!.(peps.A, Inf)
-    env, = leading_boundary(
-        CTMRGEnv(wts, peps), peps;
-        tol = ctmrg_tol, trunc = trunc_env
-    )
-    e_site = cost_function(peps, env, H) / (Nr * Nc)
+    env = CTMRGEnv(wts, peps)
+    env, = leading_boundary(env, peps; tol = ctmrg_tol, trunc = trunc_env)
+    e_site = cost_function(peps, env, ham) / (Nr * Nc)
     @info "2-site simple update energy = $e_site"
     # continue with 3-site simple update; energy should not change much
     dts = [1.0e-2, 5.0e-3]
@@ -139,14 +137,14 @@ end
     trunc = truncerror(; atol = 1.0e-10) & truncrank(2)
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
         alg = SimpleUpdate(;
-            ψ0 = peps, env0 = wts, H, dt, nstep, tol, trunc,
+            ψ0 = peps, env0 = wts, H = ham, dt, nstep, tol, trunc,
             check_interval = 1000, force_3site = true
         )
         peps, wts, = time_evolve(alg)
     end
     normalize!.(peps.A, Inf)
     env, = leading_boundary(env, peps; tol = ctmrg_tol, trunc = trunc_env)
-    e_site2 = cost_function(peps, env, H) / (Nr * Nc)
+    e_site2 = cost_function(peps, env, ham) / (Nr * Nc)
     @info "3-site simple update energy = $e_site2"
     @test e_site ≈ e_site2 atol = 5.0e-4
 end
diff --git a/test/timeevol/sitedep_truncation.jl b/test/timeevol/sitedep_truncation.jl
index 3f7a48ad..b872b0a0 100644
--- a/test/timeevol/sitedep_truncation.jl
+++ b/test/timeevol/sitedep_truncation.jl
@@ -19,7 +19,7 @@ end
 
 @testset "Simple update: bipartite 2-site" begin
     Nr, Nc = 2, 2
-    H = real(heisenberg_XYZ(InfiniteSquare(Nr, Nc); Jx = 1.0, Jy = 1.0, Jz = 1.0))
+    ham = real(heisenberg_XYZ(InfiniteSquare(Nr, Nc); Jx = 1.0, Jy = 1.0, Jz = 1.0))
     Random.seed!(100)
     peps0 = InfinitePEPS(rand, Float64, ℂ^2, ℂ^10; unitcell = (Nr, Nc))
     env0 = SUWeight(peps0)
@@ -27,7 +27,7 @@ end
     # set trunc to be compatible with bipartite structure
     bonddims = stack([[6 4; 4 6], [5 7; 7 5]]; dims = 1)
     trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims))
-    alg = SimpleUpdate(; ψ0 = peps0, env0, H, dt = 1.0e-2, nstep = 4, trunc, bipartite = true)
+    alg = SimpleUpdate(; ψ0 = peps0, env0, H = ham, dt = 1.0e-2, nstep = 4, trunc, bipartite = true)
     peps, env, = time_evolve(alg)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
@@ -44,6 +44,7 @@ end
 
 @testset "Simple update: generic 2-site and 3-site" begin
     Nr, Nc = 3, 4
+    ham = real(heisenberg_XYZ(InfiniteSquare(Nr, Nc); Jx = 1.0, Jy = 1.0, Jz = 1.0))
     Random.seed!(100)
     peps0 = InfinitePEPS(rand, Float64, ℂ^2, ℂ^10; unitcell = (Nr, Nc))
     normalize!.(peps0.A, Inf)
@@ -53,14 +54,12 @@ end
     @show bonddims
     trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims))
     # 2-site SU
-    H = real(heisenberg_XYZ(InfiniteSquare(Nr, Nc); Jx = 1.0, Jy = 1.0, Jz = 1.0))
-    alg = SimpleUpdate(; ψ0 = peps0, env0, H, dt = 1.0e-2, nstep = 4, trunc)
+    alg = SimpleUpdate(; ψ0 = peps0, env0, H = ham, dt = 1.0e-2, nstep = 4, trunc)
     peps, env, = time_evolve(alg)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
     # 3-site SU
-    H = real(j1_j2_model(InfiniteSquare(Nr, Nc); J1 = 1.0, J2 = 0.5, sublattice = false))
-    alg = SimpleUpdate(; ψ0 = peps0, env0, H, dt = 1.0e-2, nstep = 4, trunc)
+    alg = SimpleUpdate(; ψ0 = peps0, env0, H = ham, dt = 1.0e-2, nstep = 4, trunc, force_3site = true)
     peps, env, = time_evolve(alg)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims

From 08c7dc8c17c9c2f091f9b7d6648af4cc809a4708 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Thu, 6 Nov 2025 10:14:11 +0800
Subject: [PATCH 11/36] Separate SU alg, state, iterator

---
 examples/heisenberg_su/main.jl                |   8 +-
 src/PEPSKit.jl                                |   2 +-
 src/algorithms/time_evolution/evoltools.jl    |   2 +-
 src/algorithms/time_evolution/simpleupdate.jl | 195 +++++++++++-------
 .../time_evolution/simpleupdate3site.jl       |   9 +-
 src/algorithms/time_evolution/time_evolve.jl  |  45 +++-
 test/bondenv/benv_fu.jl                       |   4 +-
 test/examples/j1j2_finiteT.jl                 |  20 +-
 8 files changed, 179 insertions(+), 106 deletions(-)

diff --git a/examples/heisenberg_su/main.jl b/examples/heisenberg_su/main.jl
index b8cb0eb2..1fcee9e4 100644
--- a/examples/heisenberg_su/main.jl
+++ b/examples/heisenberg_su/main.jl
@@ -72,13 +72,9 @@ dts = [1.0e-2, 1.0e-3, 4.0e-4]
 tols = [1.0e-6, 1.0e-8, 1.0e-8]
 nstep = 10000
 trunc_peps = truncerror(; atol = 1.0e-10) & truncrank(Dbond)
-
+alg = SimpleUpdate(; trunc = trunc_peps, bipartite = true, check_interval = 500)
 for (dt, tol) in zip(dts, tols)
-    alg = SimpleUpdate(;
-        ψ0 = peps, env0 = wts, H, dt, tol, nstep,
-        trunc = trunc_peps, bipartite = true, check_interval = 500
-    )
-    global peps, wts, = time_evolve(alg)
+    global peps, wts, = time_evolve(peps, H, dt, nstep, alg, wts; tol)
 end
 
 md"""
diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index eaa08ef5..2fb5131b 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -102,7 +102,7 @@ export fixedpoint
 export absorb_weight
 export ALSTruncation, FullEnvTruncation
 export SimpleUpdate
-export time_evolve
+export TimeEvolver, timestep, time_evolve
 
 export InfiniteSquareNetwork
 export InfinitePartitionFunction
diff --git a/src/algorithms/time_evolution/evoltools.jl b/src/algorithms/time_evolution/evoltools.jl
index d00555a0..d1f66bb5 100644
--- a/src/algorithms/time_evolution/evoltools.jl
+++ b/src/algorithms/time_evolution/evoltools.jl
@@ -1,6 +1,6 @@
 const InfiniteState = Union{InfinitePEPS, InfinitePEPO}
 
-function _process_timeevol_args(
+function _get_dt(
         state::InfiniteState, dt::Float64, imaginary_time::Bool
     )
     dt′ = (state isa InfinitePEPS) ? dt : (dt / 2)
diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index be56d442..51487de1 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -8,16 +8,6 @@ Algorithm struct for simple update (SU) of InfinitePEPS or InfinitePEPO.
 $(TYPEDFIELDS)
 """
 @kwdef struct SimpleUpdate <: TimeEvolution
-    # Initial state (InfinitePEPS or InfinitePEPO)
-    ψ0::InfiniteState
-    # Hamiltonian
-    H::LocalOperator
-    # Initial Bond weights
-    env0::SUWeight
-    # Trotter time step
-    dt::Float64
-    # number of iteration steps
-    nstep::Int
     # Truncation scheme after applying Trotter gates
     trunc::TruncationStrategy
     # Switch for imaginary or real time
@@ -29,12 +19,64 @@ $(TYPEDFIELDS)
     # ---- only applicable to ground state search ----
     # assume bipartite unit cell structure
     bipartite::Bool = false
-    # converged change of weight between two steps
-    tol::Float64 = 0.0
     # ---- only applicable to PEPO evolution ----
     gate_bothsides::Bool = false # when false, purified approach is assumed
 end
 
+# internal state of simple update algorithm
+struct SUState{S}
+    # number of performed iterations
+    iter::Int
+    # evolved time
+    t::Float64
+    # measure of difference (of SUWeight, energy, etc.) from last iteration
+    diff::Float64
+    # PEPS/PEPO
+    ψ::S
+    # SUWeight environment
+    env::SUWeight
+end
+
+"""
+    TimeEvolver(
+        ψ₀::InfiniteState, H::LocalOperator, dt::Float64, nstep::Int, 
+        alg::SimpleUpdate, env₀::SUWeight; 
+        tol::Float64 = 0.0, t₀::Float64 = 0.0
+    )
+
+Initialize a TimeEvolver with Hamiltonian `H` and simple update `alg`, 
+starting from the initial state `ψ₀` and SUWeight environment `env₀`.
+The initial time is specified by `t₀`.
+Convergence check is enabled by setting `tol > 0`; otherwise it is disabled.
+"""
+function TimeEvolver(
+        ψ₀::InfiniteState, H::LocalOperator, dt::Float64, nstep::Int,
+        alg::SimpleUpdate, env₀::SUWeight;
+        tol::Float64 = 0.0, t₀::Float64 = 0.0
+    )
+    # sanity checks
+    _timeevol_sanity_check(ψ₀, H, tol, alg)
+    # create Trotter gates
+    nnonly = is_nearest_neighbour(H)
+    use_3site = alg.force_3site || !nnonly
+    if alg.bipartite
+        @assert !use_3site "3-site simple update is incompatible with bipartite lattice."
+    end
+    dt′ = _get_dt(ψ₀, dt, alg.imaginary_time)
+    gate = if use_3site
+        [
+            _get_gatempos_se(H, dt′),
+            _get_gatempos_se(rotl90(H), dt′),
+            _get_gatempos_se(rot180(H), dt′),
+            _get_gatempos_se(rotr90(H), dt′),
+        ]
+    else
+        get_expham(H, dt′)
+    end
+    state = SUState(0, t₀, NaN, ψ₀, env₀)
+    return TimeEvolver(alg, dt, nstep, H, gate, tol, state)
+end
+
 #=
 Simple update of the x-bond between `[r,c]` and `[r,c+1]`.
 
@@ -48,8 +90,7 @@ function _su_xbond!(
         state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
         row::Int, col::Int, trunc::TruncationStrategy; gate_ax::Int = 1
     ) where {T <: Number, S <: ElementarySpace}
-    Nr, Nc = size(state)[1:2]
-    @assert 1 <= row <= Nr && 1 <= col <= Nc
+    Nr, Nc, = size(state)
     cp1 = _next(col, Nc)
     # absorb environment weights
     A, B = state.A[row, col], state.A[row, cp1]
@@ -87,8 +128,7 @@ function _su_ybond!(
         state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
         row::Int, col::Int, trunc::TruncationStrategy; gate_ax::Int = 1
     ) where {T <: Number, S <: ElementarySpace}
-    Nr, Nc = size(state)[1:2]
-    @assert 1 <= row <= Nr && 1 <= col <= Nc
+    Nr, Nc, = size(state)
     rm1 = _prev(row, Nr)
     # absorb environment weights
     A, B = state.A[row, col], state.A[rm1, col]
@@ -115,19 +155,15 @@ end
 function su_iter(
         state::InfiniteState, gate::LocalOperator, alg::SimpleUpdate, env::SUWeight
     )
-    @assert size(gate.lattice) == size(state)[1:2]
-    Nr, Nc = size(state)[1:2]
-    alg.bipartite && (@assert Nr == Nc == 2)
-    (Nr >= 2 && Nc >= 2) || throw(
-        ArgumentError("`state` unit cell size for simple update should be no smaller than (2, 2)."),
-    )
-    state2, env2 = deepcopy(state), deepcopy(env)
+    Nr, Nc, = size(state)
+    state2, env2, ϵ = deepcopy(state), deepcopy(env), 0.0
     gate_axs = alg.gate_bothsides ? (1:2) : (1:1)
     for r in 1:Nr, c in 1:Nc
         term = get_gateterm(gate, (CartesianIndex(r, c), CartesianIndex(r, c + 1)))
         trunc = truncation_strategy(alg.trunc, 1, r, c)
         for gate_ax in gate_axs
-            _su_xbond!(state2, term, env2, r, c, trunc; gate_ax)
+            ϵ′ = _su_xbond!(state2, term, env2, r, c, trunc; gate_ax)
+            ϵ = maximum(ϵ, ϵ′)
         end
         if alg.bipartite
             rp1, cp1 = _next(r, Nr), _next(c, Nc)
@@ -138,7 +174,8 @@ function su_iter(
         term = get_gateterm(gate, (CartesianIndex(r, c), CartesianIndex(r - 1, c)))
         trunc = truncation_strategy(alg.trunc, 2, r, c)
         for gate_ax in gate_axs
-            _su_ybond!(state2, term, env2, r, c, trunc; gate_ax)
+            ϵ′ = _su_ybond!(state2, term, env2, r, c, trunc; gate_ax)
+            ϵ = maximum(ϵ, ϵ′)
         end
         if alg.bipartite
             rm1, cm1 = _prev(r, Nr), _prev(c, Nc)
@@ -147,61 +184,75 @@ function su_iter(
             env2.data[2, rm1, cm1] = deepcopy(env2.data[2, r, c])
         end
     end
-    wtdiff = compare_weights(env2, env)
-    return state2, env2, (; wtdiff)
+    return state2, env2, ϵ
 end
 
-# Initial iteration
-function Base.iterate(alg::SimpleUpdate)
-    # sanity checks
-    nnonly = is_nearest_neighbour(alg.H)
-    use_3site = alg.force_3site || !nnonly
-    @assert alg.tol >= 0
-    @assert !(alg.bipartite && alg.ψ0 isa InfinitePEPO) "Evolution of PEPO with bipartite structure is not implemented."
-    @assert !(alg.bipartite && use_3site) "3-site simple update is incompatible with bipartite lattice."
-    if alg.ψ0 isa InfinitePEPS
-        @assert !(alg.gate_bothsides) "alg.gate_bothsides = true is incompatible with PEPS evolution."
-    end
-    # construct Trotter 2-site gates or 3-site gate-MPOs
-    dt′ = _process_timeevol_args(alg.ψ0, alg.dt, alg.imaginary_time)
-    gate = if use_3site
-        [
-            _get_gatempos_se(alg.H, dt′),
-            _get_gatempos_se(rotl90(alg.H), dt′),
-            _get_gatempos_se(rot180(alg.H), dt′),
-            _get_gatempos_se(rotr90(alg.H), dt′),
-        ]
+function Base.iterate(it::TimeEvolver{<:SimpleUpdate}, state = it.state)
+    alg = it.alg
+    iter, t, diff = state.iter, state.t, state.diff
+    check_convergence = (it.tol > 0)
+    if check_convergence
+        if diff < it.tol
+            @info "SU: bond weights have converged."
+            return nothing
+        elseif iter == it.nstep
+            @warn "SU: bond weights have not converged."
+            return nothing
+        end
     else
-        get_expham(alg.H, dt′)
+        (iter == it.nstep) && return nothing
     end
+    # perform one iteration and record time
     time0 = time()
-    ψ, env, info = su_iter(alg.ψ0, gate, alg, alg.env0)
+    ψ, env, ϵ = su_iter(state.ψ, it.gate, it.alg, state.env)
+    diff = compare_weights(env, state.env)
     elapsed_time = time() - time0
-    converged = false
-    return (ψ, env, info), (ψ, env, info, gate, converged, 1, elapsed_time)
-end
-
-# subsequent iterations
-function Base.iterate(alg::SimpleUpdate, state)
-    ψ0, env0, info, gate, converged, it, elapsed_time = state
-    check_convergence = (alg.tol > 0)
-    if (it % alg.check_interval == 0) || (it == 1) || (it == alg.nstep) || converged
-        @info "Space of x-weight at [1, 1] = $(space(env0[1, 1, 1], 1))"
+    # update internal state
+    iter += 1
+    t += it.dt
+    it.state = SUState(iter, t, diff, ψ, env)
+    # output information
+    converged = check_convergence ? (diff < it.tol) : false
+    showinfo = (iter % alg.check_interval == 0) || (iter == 1) || (iter == it.nstep) || converged
+    if showinfo
+        @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))"
         @info @sprintf(
             "SU iter %-7d: dt = %.0e, |Δλ| = %.3e. Time = %.3f s/it",
-            it, alg.dt, info.wtdiff, elapsed_time
+            iter, it.dt, diff, elapsed_time
         )
-        if (it == alg.nstep)
-            check_convergence && (@warn "SU: bond weights have not converged.")
-            return nothing
-        elseif converged
-            @info "SU: bond weights have converged."
-            return nothing
-        end
     end
-    time0 = time()
-    ψ, env, info = su_iter(ψ0, gate, alg, env0)
-    elapsed_time = time() - time0
-    converged = check_convergence ? (info.wtdiff < alg.tol) : false
-    return (ψ, env, info), (ψ, env, info, gate, converged, it + 1, elapsed_time)
+    info = (; t, ϵ)
+    return (ψ, env, info), it.state
+end
+
+# evolve one step from given state
+function MPSKit.timestep(
+        it::TimeEvolver{<:SimpleUpdate}, ψ::InfiniteState, env::SUWeight;
+        iter::Int = it.state.iter, t::Float64 = it.state.t
+    )
+    _timeevol_sanity_check(ψ, it.H, it.tol, it.alg)
+    state = SUState(iter, t, NaN, ψ, env)
+    return first(iterate(it, state))
+end
+
+# evolve till the end
+function MPSKit.time_evolve(it::TimeEvolver{<:SimpleUpdate})
+    time_start = time()
+    result = nothing
+    for state in it
+        result = state
+    end
+    time_end = time()
+    @info @sprintf("Simple update finished. Total time elasped: %.2f s", time_end - time_start)
+    return result
+end
+
+# ordinary mode
+function MPSKit.time_evolve(
+        ψ₀::InfiniteState, H::LocalOperator, dt::Float64, nstep::Int,
+        alg::SimpleUpdate, env₀::SUWeight;
+        tol::Float64 = 0.0, t₀::Float64 = 0.0
+    )
+    it = TimeEvolver(ψ₀, H, dt, nstep, alg, env₀; tol, t₀)
+    return time_evolve(it)
 end
diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl
index 7dcbdd08..c75663b9 100644
--- a/src/algorithms/time_evolution/simpleupdate3site.jl
+++ b/src/algorithms/time_evolution/simpleupdate3site.jl
@@ -500,7 +500,7 @@ function _su3site_se!(
         M = absorb_weight(M, env, coord[1], coord[2], openaxs; inv = true)
         state.A[CartesianIndex(coord)] = normalize(M, Inf)
     end
-    return ϵs
+    return maximum(ϵs)
 end
 
 function su_iter(
@@ -515,7 +515,7 @@ function su_iter(
             "iPEPS unit cell size for simple update should be no smaller than (2, 2)."
         ),
     )
-    state2, env2 = deepcopy(state), deepcopy(env)
+    state2, env2, ϵ = deepcopy(state), deepcopy(env), 0.0
     trunc = alg.trunc
     for i in 1:4
         Nr, Nc = size(state2)[1:2]
@@ -525,11 +525,10 @@ function su_iter(
                 truncation_strategy(trunc, 1, r, c)
                 truncation_strategy(trunc, 2, r, _next(c, Nc))
             ]
-            _su3site_se!(state2, gs, env2, r, c, truncs; alg.gate_bothsides)
+            ϵ = _su3site_se!(state2, gs, env2, r, c, truncs; alg.gate_bothsides)
         end
         state2, env2 = rotl90(state2), rotl90(env2)
         trunc = rotl90(trunc)
     end
-    wtdiff = compare_weights(env2, env)
-    return state2, env2, (; wtdiff)
+    return state2, env2, ϵ
 end
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index 6ef53339..e5bd0ed9 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -1,12 +1,41 @@
 abstract type TimeEvolution end
 
-function MPSKit.time_evolve(alg::Alg) where {Alg <: TimeEvolution}
-    time_start = time()
-    result = nothing
-    for state in alg
-        result = state
+mutable struct TimeEvolver{TE <: TimeEvolution, G, S}
+    # Time evolution algorithm
+    alg::TE
+    # Trotter time step
+    dt::Float64
+    # Maximal iteration steps
+    nstep::Int
+    # Hamiltonian
+    H::LocalOperator
+    # Trotter gates
+    gate::G
+    # Convergence tolerance (change of weight or energy from last iteration)
+    tol::Float64
+    # PEPS/PEPO (and environment)
+    state::S
+end
+
+Base.iterate(it::TimeEvolver) = iterate(it, it.state)
+
+function _timeevol_sanity_check(
+        ψ₀::InfiniteState, H::LocalOperator, tol::Float64, alg::A
+    ) where {A <: TimeEvolution}
+    Nr, Nc, = size(ψ₀)
+    @assert (Nr >= 2 && Nc >= 2) "Unit cell size for simple update should be no smaller than (2, 2)."
+    @assert physicalspace(H) == physicalspace(ψ₀) "Physical space mismatch between `ψ₀` and `H`."
+    @assert tol >= 0
+    if tol > 0
+        @assert alg.imaginary_time "`tol` should be 0 for real time evolution."
+        @assert ψ₀ isa InfinitePEPS "`tol` should be 0 for time evolution of InfinitePEPO."
+    end
+    if hasfield(typeof(alg), :gate_bothsides) && alg.gate_bothsides
+        @assert ψ₀ isa InfinitePEPO "alg.gate_bothsides = true is only compatible with PEPO."
+    end
+    if hasfield(typeof(alg), :bipartite) && alg.bipartite
+        @assert Nr == Nc == 2 "`bipartite = true` requires 2 x 2 unit cell size."
+        @assert ψ₀ isa InfinitePEPS "Evolution of PEPO with bipartite structure is not implemented."
     end
-    time_end = time()
-    @info @sprintf("Time evolution finished. Total time elasped: %.2f s", time_end - time_start)
-    return result
+    return nothing
 end
diff --git a/test/bondenv/benv_fu.jl b/test/bondenv/benv_fu.jl
index 181616b1..8d74f96a 100644
--- a/test/bondenv/benv_fu.jl
+++ b/test/bondenv/benv_fu.jl
@@ -14,10 +14,10 @@ function get_hubbard_state(t::Float64 = 1.0, U::Float64 = 8.0)
     peps = InfinitePEPS(rand, ComplexF64, Vphy, Vphy; unitcell = (Nr, Nc))
     wts = SUWeight(peps)
     alg = SimpleUpdate(;
-        ψ0 = peps, env0 = wts, H, dt = 1.0e-2, nstep = 10000, tol = 1.0e-8,
         trunc = truncerror(; atol = 1.0e-10) & truncrank(4), check_interval = 2000
     )
-    peps, = time_evolve(alg)
+    evolver = TimeEvolver(peps, H, 1.0e-2, 10000, alg, wts; tol = 1.0e-8)
+    peps, = time_evolve(evolver)
     normalize!.(peps.A, Inf)
     return peps
 end
diff --git a/test/examples/j1j2_finiteT.jl b/test/examples/j1j2_finiteT.jl
index 9a3828a8..8ffa3836 100644
--- a/test/examples/j1j2_finiteT.jl
+++ b/test/examples/j1j2_finiteT.jl
@@ -32,22 +32,19 @@ check_interval = 100
 
 # PEPO approach
 dt, nstep = 1.0e-3, 600
-alg = SimpleUpdate(;
-    ψ0 = pepo0, env0 = wts0, H = ham, dt, nstep, trunc = trunc_pepo,
-    gate_bothsides = true, check_interval
-)
-pepo, wts, = time_evolve(alg)
+alg = SimpleUpdate(; trunc = trunc_pepo, gate_bothsides = true, check_interval)
+evolver = TimeEvolver(pepo0, ham, dt, nstep, alg, wts0)
+pepo, wts, info = time_evolve(evolver)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 energy = expectation_value(pepo, ham, env) / (Nr * Nc)
 @info "β = $(dt * nstep): tr(ρH) = $(energy)"
+@test dt * nstep ≈ info.t
 @test energy ≈ bm[2] atol = 5.0e-3
 
 # PEPS (purified PEPO) approach
-alg = SimpleUpdate(;
-    ψ0 = pepo0, env0 = wts0, H = ham, dt, nstep, trunc = trunc_pepo,
-    gate_bothsides = false, check_interval
-)
-pepo, wts, = time_evolve(alg)
+alg = SimpleUpdate(; trunc = trunc_pepo, gate_bothsides = false, check_interval)
+evolver = TimeEvolver(pepo0, ham, dt, nstep, alg, wts0)
+pepo, wts, = time_evolve(evolver)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 energy = expectation_value(pepo, ham, env) / (Nr * Nc)
 @info "β = $(dt * nstep) / 2: tr(ρH) = $(energy)"
@@ -55,5 +52,6 @@ energy = expectation_value(pepo, ham, env) / (Nr * Nc)
 
 env = converge_env(InfinitePEPS(pepo), 16)
 energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc)
-@info "β = $(dt * nstep): ⟨ρ|ham|ρ⟩ = $(energy)"
+@info "β = $(dt * nstep): ⟨ρ|H|ρ⟩ = $(energy)"
+@test dt * nstep ≈ info.t
 @test energy ≈ bm[2] atol = 5.0e-3

From d8c2653908c9b465987234045e4c8661db2d9966 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Thu, 6 Nov 2025 11:10:43 +0800
Subject: [PATCH 12/36] Remove `H` from TimeEvolver

---
 src/algorithms/time_evolution/simpleupdate.jl | 11 ++++++-----
 src/algorithms/time_evolution/time_evolve.jl  |  8 +++-----
 2 files changed, 9 insertions(+), 10 deletions(-)

diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index 51487de1..5cebb080 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -55,7 +55,7 @@ function TimeEvolver(
         tol::Float64 = 0.0, t₀::Float64 = 0.0
     )
     # sanity checks
-    _timeevol_sanity_check(ψ₀, H, tol, alg)
+    _timeevol_sanity_check(ψ₀, physicalspace(H), tol, alg)
     # create Trotter gates
     nnonly = is_nearest_neighbour(H)
     use_3site = alg.force_3site || !nnonly
@@ -74,7 +74,7 @@ function TimeEvolver(
         get_expham(H, dt′)
     end
     state = SUState(0, t₀, NaN, ψ₀, env₀)
-    return TimeEvolver(alg, dt, nstep, H, gate, tol, state)
+    return TimeEvolver(alg, dt, nstep, gate, tol, state)
 end
 
 #=
@@ -163,7 +163,7 @@ function su_iter(
         trunc = truncation_strategy(alg.trunc, 1, r, c)
         for gate_ax in gate_axs
             ϵ′ = _su_xbond!(state2, term, env2, r, c, trunc; gate_ax)
-            ϵ = maximum(ϵ, ϵ′)
+            ϵ = max(ϵ, ϵ′)
         end
         if alg.bipartite
             rp1, cp1 = _next(r, Nr), _next(c, Nc)
@@ -175,7 +175,7 @@ function su_iter(
         trunc = truncation_strategy(alg.trunc, 2, r, c)
         for gate_ax in gate_axs
             ϵ′ = _su_ybond!(state2, term, env2, r, c, trunc; gate_ax)
-            ϵ = maximum(ϵ, ϵ′)
+            ϵ = max(ϵ, ϵ′)
         end
         if alg.bipartite
             rm1, cm1 = _prev(r, Nr), _prev(c, Nc)
@@ -230,7 +230,8 @@ function MPSKit.timestep(
         it::TimeEvolver{<:SimpleUpdate}, ψ::InfiniteState, env::SUWeight;
         iter::Int = it.state.iter, t::Float64 = it.state.t
     )
-    _timeevol_sanity_check(ψ, it.H, it.tol, it.alg)
+    Pspaces = physicalspace(it.state.ψ)
+    _timeevol_sanity_check(ψ, Pspaces, it.tol, it.alg)
     state = SUState(iter, t, NaN, ψ, env)
     return first(iterate(it, state))
 end
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index e5bd0ed9..af51a714 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -7,8 +7,6 @@ mutable struct TimeEvolver{TE <: TimeEvolution, G, S}
     dt::Float64
     # Maximal iteration steps
     nstep::Int
-    # Hamiltonian
-    H::LocalOperator
     # Trotter gates
     gate::G
     # Convergence tolerance (change of weight or energy from last iteration)
@@ -20,11 +18,11 @@ end
 Base.iterate(it::TimeEvolver) = iterate(it, it.state)
 
 function _timeevol_sanity_check(
-        ψ₀::InfiniteState, H::LocalOperator, tol::Float64, alg::A
-    ) where {A <: TimeEvolution}
+        ψ₀::InfiniteState, Pspaces::M, tol::Float64, alg::A
+    ) where {A <: TimeEvolution, M <: AbstractMatrix{<:ElementarySpace}}
     Nr, Nc, = size(ψ₀)
     @assert (Nr >= 2 && Nc >= 2) "Unit cell size for simple update should be no smaller than (2, 2)."
-    @assert physicalspace(H) == physicalspace(ψ₀) "Physical space mismatch between `ψ₀` and `H`."
+    @assert Pspaces == physicalspace(ψ₀) "Physical spaces of `ψ₀` do not match `Pspaces`."
     @assert tol >= 0
     if tol > 0
         @assert alg.imaginary_time "`tol` should be 0 for real time evolution."

From f589d9d4fb47992db55534b19d704ce9b53c09f5 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Thu, 6 Nov 2025 11:44:28 +0800
Subject: [PATCH 13/36] Add test for `timestep`

---
 src/algorithms/time_evolution/simpleupdate.jl |  8 ++++-
 test/runtests.jl                              |  3 ++
 test/timeevol/timestep.jl                     | 33 +++++++++++++++++++
 3 files changed, 43 insertions(+), 1 deletion(-)
 create mode 100644 test/timeevol/timestep.jl

diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index 5cebb080..cd601089 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -233,7 +233,13 @@ function MPSKit.timestep(
     Pspaces = physicalspace(it.state.ψ)
     _timeevol_sanity_check(ψ, Pspaces, it.tol, it.alg)
     state = SUState(iter, t, NaN, ψ, env)
-    return first(iterate(it, state))
+    result = iterate(it, state)
+    if result === nothing
+        @warn "TimeEvolver `it` has already reached the end."
+        return ψ, env, (; t, ϵ = NaN)
+    else
+        return first(result)
+    end
 end
 
 # evolve till the end
diff --git a/test/runtests.jl b/test/runtests.jl
index f1800794..38d4fdd4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -61,6 +61,9 @@ end
         end
     end
     if GROUP == "ALL" || GROUP == "TIMEEVOL"
+        @time @safetestset "`timestep` function" begin
+            include("timeevol/timestep.jl")
+        end
         @time @safetestset "Cluster truncation with projectors" begin
             include("timeevol/cluster_projectors.jl")
         end
diff --git a/test/timeevol/timestep.jl b/test/timeevol/timestep.jl
new file mode 100644
index 00000000..ede88db6
--- /dev/null
+++ b/test/timeevol/timestep.jl
@@ -0,0 +1,33 @@
+using Test
+using Random
+using TensorKit
+using PEPSKit
+
+function test_timestep(ψ0, env0, H)
+    trunc = truncerror(; atol = 1.0e-10) & truncrank(4)
+    alg = SimpleUpdate(; trunc)
+    dt, nstep = 1.0e-2, 50
+    # manual timestep
+    evolver = TimeEvolver(ψ0, H, dt, nstep, alg, env0)
+    ψ1, env1, info1 = deepcopy(ψ0), deepcopy(env0), nothing
+    for iter in 0:(nstep - 1)
+        ψ1, env1, info1 = timestep(evolver, ψ1, env1; iter)
+    end
+    @info info1
+    # time_evolve
+    ψ2, env2, info2 = time_evolve(ψ0, H, dt, nstep, alg, env0)
+    @info info2
+    # results should be *exactly* the same
+    @test ψ1 == ψ2
+    @test env1 == env2
+    @test info1 == info2
+    return nothing
+end
+
+Nr, Nc = 2, 2
+H = real(heisenberg_XYZ(ComplexF64, Trivial, InfiniteSquare(Nr, Nc); Jx = 1, Jy = 1, Jz = 1))
+Pspace, Vspace = ℂ^2, ℂ^4
+ψ0 = InfinitePEPS(rand, Float64, Pspace, Vspace; unitcell = (Nr, Nc))
+env0 = SUWeight(ψ0)
+
+test_timestep(ψ0, env0, H)

From 2b8fa1fdd26d028849deefbad77de94cd9f63eba Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Thu, 6 Nov 2025 12:04:18 +0800
Subject: [PATCH 14/36] Add docstrings

---
 src/algorithms/time_evolution/simpleupdate.jl | 44 ++++++++++++++++---
 src/algorithms/time_evolution/time_evolve.jl  | 14 ++++++
 2 files changed, 53 insertions(+), 5 deletions(-)

diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index cd601089..8877ce85 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -46,8 +46,10 @@ end
 
 Initialize a TimeEvolver with Hamiltonian `H` and simple update `alg`, 
 starting from the initial state `ψ₀` and SUWeight environment `env₀`.
-The initial time is specified by `t₀`.
-Convergence check is enabled by setting `tol > 0`; otherwise it is disabled.
+
+- The initial (real or imaginary) time is specified by `t₀`.
+- Setting `tol > 0` enables convergence check (for imaginary time evolution of InfinitePEPS only).
+    For other usages it should not be changed.
 """
 function TimeEvolver(
         ψ₀::InfiniteState, H::LocalOperator, dt::Float64, nstep::Int,
@@ -225,7 +227,18 @@ function Base.iterate(it::TimeEvolver{<:SimpleUpdate}, state = it.state)
     return (ψ, env, info), it.state
 end
 
-# evolve one step from given state
+"""
+    timestep(
+        it::TimeEvolver{<:SimpleUpdate}, ψ::InfiniteState, env::SUWeight;
+        iter::Int = it.state.iter, t::Float64 = it.state.t
+    ) -> (ψ, env, info)
+
+Given the TimeEvolver iterator `it`, perform one step of time evolution
+on the input state `ψ` and its environment `env`.
+
+- Using `iter` and `t` to reset the current iteration number and evolved time
+    respectively of the TimeEvolver `it`.
+"""
 function MPSKit.timestep(
         it::TimeEvolver{<:SimpleUpdate}, ψ::InfiniteState, env::SUWeight;
         iter::Int = it.state.iter, t::Float64 = it.state.t
@@ -242,7 +255,12 @@ function MPSKit.timestep(
     end
 end
 
-# evolve till the end
+"""
+    time_evolve(it::TimeEvolver{<:SimpleUpdate})
+
+Perform time evolution until the set number of iterations or convergence
+directly using the specified TimeEvolver iterator.
+"""
 function MPSKit.time_evolve(it::TimeEvolver{<:SimpleUpdate})
     time_start = time()
     result = nothing
@@ -254,7 +272,23 @@ function MPSKit.time_evolve(it::TimeEvolver{<:SimpleUpdate})
     return result
 end
 
-# ordinary mode
+"""
+    time_evolve(
+        ψ₀::Union{InfinitePEPS, InfinitePEPO}, H::LocalOperator, 
+        dt::Float64, nstep::Int, alg::SimpleUpdate, env₀::SUWeight;
+        tol::Float64 = 0.0, t₀::Float64 = 0.0
+    ) -> (ψ, env, info)
+
+Perform time evolution on the initial state `ψ₀` and initial environment `env₀`
+with Hamiltonian `H`, using SimpleUpdate algorithm `alg`, time step `dt` for 
+`nstep` number of steps. 
+
+- Setting `tol > 0` enables convergence check (for imaginary time evolution of InfinitePEPS only).
+    For other usages it should not be changed.
+- Using `t₀` to specify the initial (real or imaginary) time of `ψ₀`.
+- `info` is a NamedTuple containing information of the evolution, 
+    including the time evolved since `ψ₀`.
+"""
 function MPSKit.time_evolve(
         ψ₀::InfiniteState, H::LocalOperator, dt::Float64, nstep::Int,
         alg::SimpleUpdate, env₀::SUWeight;
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index af51a714..cb4ae3ce 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -1,5 +1,19 @@
+"""
+$(TYPEDEF)
+
+Abstract super type for time evolution algorithms of InfinitePEPS or InfinitePEPO.
+"""
 abstract type TimeEvolution end
 
+"""
+    mutable struct TimeEvolver{TE <: TimeEvolution, G, S}
+
+Iterator for Trotter-based time evolution of InfinitePEPS or InfinitePEPO.
+
+## Fields
+
+$(TYPEDFIELDS)
+"""
 mutable struct TimeEvolver{TE <: TimeEvolution, G, S}
     # Time evolution algorithm
     alg::TE

From 0844b6081bb98ecf393ad881c8186bcb4938ded9 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Thu, 6 Nov 2025 12:40:48 +0800
Subject: [PATCH 15/36] Update SU examples

---
 examples/hubbard_su/main.jl                   | 16 +++++++---------
 examples/j1j2_su/main.jl                      | 14 +++++---------
 src/algorithms/time_evolution/simpleupdate.jl |  4 ++--
 3 files changed, 14 insertions(+), 20 deletions(-)

diff --git a/examples/hubbard_su/main.jl b/examples/hubbard_su/main.jl
index 68c94cc2..bfb5350a 100644
--- a/examples/hubbard_su/main.jl
+++ b/examples/hubbard_su/main.jl
@@ -66,26 +66,24 @@ maxiter = 20000
 
 for (dt, tol, Dbond) in zip(dts, tols, Ds)
     trunc = truncerror(; atol = 1.0e-10) & truncrank(Dbond)
-    alg = SimpleUpdate(dt, tol, maxiter, trunc)
-    global peps, wts, = simpleupdate(
-        peps, H, alg, wts; bipartite = false, check_interval = 2000
-    )
+    alg = SimpleUpdate(; trunc, bipartite = false, check_interval = 2000)
+    global peps, wts, = time_evolve(peps, H, dt, maxiter, alg, wts; tol)
 end
 
 md"""
 ## Computing the ground-state energy
 
 In order to compute the energy expectation value with evolved PEPS, we need to converge a
-CTMRG environment on it. We first converge an environment with a small enviroment dimension
-and then use that to initialize another run with bigger environment dimension. We'll use
-`trunc=truncrank(χ)` for that such that the dimension is increased during the second CTMRG
-run:
+CTMRG environment on it. We first converge an environment with a small enviroment dimension,
+which is initialized using the simple update bond weights. Next we use it to initialize 
+another run with bigger environment dimension. The dynamic adjustion of environment dimension
+is achieved by using `trunc=truncrank(χ)` in the CTMRG runs:
 """
 
 χenv₀, χenv = 6, 16
 env_space = Vect[fℤ₂](0 => χenv₀ / 2, 1 => χenv₀ / 2)
 normalize!.(peps.A, Inf)
-env = CTMRGEnv(rand, Float64, peps, env_space)
+env = CTMRGEnv(wts, peps)
 for χ in [χenv₀, χenv]
     global env, = leading_boundary(
         env, peps; alg = :sequential, tol = 1.0e-8, maxiter = 50, trunc = truncrank(χ)
diff --git a/examples/j1j2_su/main.jl b/examples/j1j2_su/main.jl
index d50e256e..6b3ec575 100644
--- a/examples/j1j2_su/main.jl
+++ b/examples/j1j2_su/main.jl
@@ -51,14 +51,13 @@ on the previously evolved PEPS:
 dt, tol, nstep = 1.0e-2, 1.0e-8, 30000
 check_interval = 4000
 trunc_peps = truncerror(; atol = 1.0e-10) & truncrank(Dbond)
+alg = SimpleUpdate(; trunc = trunc_peps, check_interval)
 for J2 in 0.1:0.1:0.5
-    H = real( ## convert Hamiltonian `LocalOperator` to real floats
+    ## convert Hamiltonian `LocalOperator` to real floats
+    H = real(
         j1_j2_model(ComplexF64, symm, InfiniteSquare(Nr, Nc); J1, J2, sublattice = false),
     )
-    alg = SimpleUpdate(;
-        ψ0 = peps, env0 = wts, H, dt, tol, nstep, trunc = trunc_peps, check_interval
-    )
-    global peps, wts, = time_evolve(alg)
+    global peps, wts, = time_evolve(peps, H, dt, nstep, alg, wts; tol)
 end
 
 md"""
@@ -71,10 +70,7 @@ tols = [1.0e-9, 1.0e-9]
 J2 = 0.5
 H = real(j1_j2_model(ComplexF64, symm, InfiniteSquare(Nr, Nc); J1, J2, sublattice = false))
 for (dt, tol) in zip(dts, tols)
-    alg = SimpleUpdate(;
-        ψ0 = peps, env0 = wts, H, dt, tol, nstep, trunc = trunc_peps, check_interval
-    )
-    global peps, wts, = time_evolve(alg)
+    global peps, wts, = time_evolve(peps, H, dt, nstep, alg, wts; tol)
 end
 
 md"""
diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index 8877ce85..da0b352a 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -258,8 +258,8 @@ end
 """
     time_evolve(it::TimeEvolver{<:SimpleUpdate})
 
-Perform time evolution until the set number of iterations or convergence
-directly using the specified TimeEvolver iterator.
+Perform time evolution until convergence or the set number of iterations
+using the specified TimeEvolver iterator `it` directly.
 """
 function MPSKit.time_evolve(it::TimeEvolver{<:SimpleUpdate})
     time_start = time()

From 9e1f8252a6b5de59908bbc2af00c2665ec43fc60 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Thu, 6 Nov 2025 14:56:22 +0800
Subject: [PATCH 16/36] Update tests

---
 test/examples/heisenberg.jl         |  4 ++--
 test/examples/tf_ising_finiteT.jl   | 23 ++++++++---------------
 test/timeevol/cluster_projectors.jl | 14 ++++----------
 test/timeevol/sitedep_truncation.jl | 12 ++++++------
 4 files changed, 20 insertions(+), 33 deletions(-)

diff --git a/test/examples/heisenberg.jl b/test/examples/heisenberg.jl
index d4224af9..df526e9f 100644
--- a/test/examples/heisenberg.jl
+++ b/test/examples/heisenberg.jl
@@ -71,8 +71,8 @@ end
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
         Dbond2 = (n == 2) ? Dbond + 2 : Dbond
         trunc = truncerror(; atol = 1.0e-10) & truncrank(Dbond2)
-        alg = SimpleUpdate(; ψ0 = peps, env0 = wts, H = ham, dt, nstep, trunc, bipartite = false)
-        peps, wts, = time_evolve(alg)
+        alg = SimpleUpdate(; trunc, bipartite = false)
+        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts; tol)
     end
 
     # measure physical quantities with CTMRG
diff --git a/test/examples/tf_ising_finiteT.jl b/test/examples/tf_ising_finiteT.jl
index f3217b90..625a209e 100644
--- a/test/examples/tf_ising_finiteT.jl
+++ b/test/examples/tf_ising_finiteT.jl
@@ -60,34 +60,27 @@ dt, nstep = 1.0e-3, 400
 # when g = 2, β = 0.4 and 2β = 0.8 belong to two phases (without and with nonzero σᶻ)
 
 # PEPO approach: results at β, or T = 2.5
-alg = SimpleUpdate(;
-    ψ0 = pepo0, env0 = wts0, H = ham, dt, nstep,
-    trunc = trunc_pepo, gate_bothsides = true
-)
-pepo, wts, = time_evolve(alg)
+alg = SimpleUpdate(; trunc = trunc_pepo, gate_bothsides = true)
+pepo, wts, info = time_evolve(pepo0, ham, dt, nstep, alg, wts0)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 result_β = measure_mag(pepo, env)
 @info "Magnetization at T = $(1 / β)" result_β
+@test β ≈ info.t
 @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2)
 
 # continue to get results at 2β, or T = 1.25
-alg = SimpleUpdate(;
-    ψ0 = pepo, env0 = wts, H = ham, dt, nstep,
-    trunc = trunc_pepo, gate_bothsides = true
-)
-pepo, wts, = time_evolve(alg)
+pepo, wts, info = time_evolve(pepo, ham, dt, nstep, alg, wts; t₀ = β)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 result_2β = measure_mag(pepo, env)
 @info "Magnetization at T = $(1 / (2β))" result_2β
+@test 2 * β ≈ info.t
 @test isapprox(abs.(result_2β), bm_2β, rtol = 1.0e-4)
 
 # Purification approach: results at 2β, or T = 1.25
-alg = SimpleUpdate(;
-    ψ0 = pepo0, env0 = wts0, H = ham, dt, nstep = 2 * nstep,
-    trunc = trunc_pepo, gate_bothsides = false
-)
-pepo, = time_evolve(alg)
+alg = SimpleUpdate(; trunc = trunc_pepo, gate_bothsides = false)
+pepo, wts, info = time_evolve(pepo0, ham, dt, 2 * nstep, alg, wts0)
 env = converge_env(InfinitePEPS(pepo), 8)
 result_2β′ = measure_mag(pepo, env; purified = true)
 @info "Magnetization at T = $(1 / (2β)) (purification approach)" result_2β′
+@test 2 * β ≈ info.t
 @test isapprox(abs.(result_2β′), bm_2β, rtol = 1.0e-2)
diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl
index 501e6334..71f4163a 100644
--- a/test/timeevol/cluster_projectors.jl
+++ b/test/timeevol/cluster_projectors.jl
@@ -120,11 +120,8 @@ end
     nstep = 5000
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
         trunc = truncerror(; atol = 1.0e-10) & truncrank(n == 1 ? 4 : 2)
-        alg = SimpleUpdate(;
-            ψ0 = peps, env0 = wts, H = ham, dt, nstep, tol, trunc,
-            bipartite = true, check_interval = 1000
-        )
-        peps, wts, = time_evolve(alg)
+        alg = SimpleUpdate(; trunc, bipartite = true, check_interval = 1000)
+        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts; tol)
     end
     normalize!.(peps.A, Inf)
     env = CTMRGEnv(wts, peps)
@@ -135,12 +132,9 @@ end
     dts = [1.0e-2, 5.0e-3]
     tols = [1.0e-8, 1.0e-8]
     trunc = truncerror(; atol = 1.0e-10) & truncrank(2)
+    alg = SimpleUpdate(; trunc, check_interval = 1000, force_3site = true)
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
-        alg = SimpleUpdate(;
-            ψ0 = peps, env0 = wts, H = ham, dt, nstep, tol, trunc,
-            check_interval = 1000, force_3site = true
-        )
-        peps, wts, = time_evolve(alg)
+        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts; tol)
     end
     normalize!.(peps.A, Inf)
     env, = leading_boundary(env, peps; tol = ctmrg_tol, trunc = trunc_env)
diff --git a/test/timeevol/sitedep_truncation.jl b/test/timeevol/sitedep_truncation.jl
index b872b0a0..fbc86fb8 100644
--- a/test/timeevol/sitedep_truncation.jl
+++ b/test/timeevol/sitedep_truncation.jl
@@ -27,8 +27,8 @@ end
     # set trunc to be compatible with bipartite structure
     bonddims = stack([[6 4; 4 6], [5 7; 7 5]]; dims = 1)
     trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims))
-    alg = SimpleUpdate(; ψ0 = peps0, env0, H = ham, dt = 1.0e-2, nstep = 4, trunc, bipartite = true)
-    peps, env, = time_evolve(alg)
+    alg = SimpleUpdate(; trunc, bipartite = true)
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
     # check bipartite structure is preserved
@@ -54,13 +54,13 @@ end
     @show bonddims
     trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims))
     # 2-site SU
-    alg = SimpleUpdate(; ψ0 = peps0, env0, H = ham, dt = 1.0e-2, nstep = 4, trunc)
-    peps, env, = time_evolve(alg)
+    alg = SimpleUpdate(; trunc)
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
     # 3-site SU
-    alg = SimpleUpdate(; ψ0 = peps0, env0, H = ham, dt = 1.0e-2, nstep = 4, trunc, force_3site = true)
-    peps, env, = time_evolve(alg)
+    alg = SimpleUpdate(; trunc, force_3site = true)
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
 end

From 3abd2ceac73abb5de078ffe63fef89f9ea878f84 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Thu, 6 Nov 2025 14:59:21 +0800
Subject: [PATCH 17/36] Minor changes

---
 test/examples/tf_ising_finiteT.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/examples/tf_ising_finiteT.jl b/test/examples/tf_ising_finiteT.jl
index 625a209e..766fb632 100644
--- a/test/examples/tf_ising_finiteT.jl
+++ b/test/examples/tf_ising_finiteT.jl
@@ -64,7 +64,7 @@ alg = SimpleUpdate(; trunc = trunc_pepo, gate_bothsides = true)
 pepo, wts, info = time_evolve(pepo0, ham, dt, nstep, alg, wts0)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 result_β = measure_mag(pepo, env)
-@info "Magnetization at T = $(1 / β)" result_β
+@info "tr(σ(x,z)ρ) at T = $(1 / β)" result_β
 @test β ≈ info.t
 @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2)
 
@@ -72,7 +72,7 @@ result_β = measure_mag(pepo, env)
 pepo, wts, info = time_evolve(pepo, ham, dt, nstep, alg, wts; t₀ = β)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 result_2β = measure_mag(pepo, env)
-@info "Magnetization at T = $(1 / (2β))" result_2β
+@info "tr(σ(x,z)ρ) at T = $(1 / (2β))" result_2β
 @test 2 * β ≈ info.t
 @test isapprox(abs.(result_2β), bm_2β, rtol = 1.0e-4)
 
@@ -81,6 +81,6 @@ alg = SimpleUpdate(; trunc = trunc_pepo, gate_bothsides = false)
 pepo, wts, info = time_evolve(pepo0, ham, dt, 2 * nstep, alg, wts0)
 env = converge_env(InfinitePEPS(pepo), 8)
 result_2β′ = measure_mag(pepo, env; purified = true)
-@info "Magnetization at T = $(1 / (2β)) (purification approach)" result_2β′
+@info "⟨ρ|σ(x,z)|ρ⟩ at T = $(1 / (2β)) (purification approach)" result_2β′
 @test 2 * β ≈ info.t
 @test isapprox(abs.(result_2β′), bm_2β, rtol = 1.0e-2)

From b5e39da85b146394f220df95e29960024dd2ddf6 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Thu, 6 Nov 2025 18:08:52 +0800
Subject: [PATCH 18/36] Add verbosity for TimeEvolver

---
 src/algorithms/time_evolution/simpleupdate.jl | 80 ++++++++++---------
 src/algorithms/time_evolution/time_evolve.jl  |  2 +
 test/timeevol/sitedep_truncation.jl           |  6 +-
 3 files changed, 46 insertions(+), 42 deletions(-)

diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index da0b352a..cdf73fb6 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -41,7 +41,7 @@ end
     TimeEvolver(
         ψ₀::InfiniteState, H::LocalOperator, dt::Float64, nstep::Int, 
         alg::SimpleUpdate, env₀::SUWeight; 
-        tol::Float64 = 0.0, t₀::Float64 = 0.0
+        tol::Float64 = 0.0, t₀::Float64 = 0.0, verbosity = 1
     )
 
 Initialize a TimeEvolver with Hamiltonian `H` and simple update `alg`, 
@@ -54,7 +54,7 @@ starting from the initial state `ψ₀` and SUWeight environment `env₀`.
 function TimeEvolver(
         ψ₀::InfiniteState, H::LocalOperator, dt::Float64, nstep::Int,
         alg::SimpleUpdate, env₀::SUWeight;
-        tol::Float64 = 0.0, t₀::Float64 = 0.0
+        tol::Float64 = 0.0, t₀::Float64 = 0.0, verbosity::Int = 1
     )
     # sanity checks
     _timeevol_sanity_check(ψ₀, physicalspace(H), tol, alg)
@@ -76,7 +76,7 @@ function TimeEvolver(
         get_expham(H, dt′)
     end
     state = SUState(0, t₀, NaN, ψ₀, env₀)
-    return TimeEvolver(alg, dt, nstep, gate, tol, state)
+    return TimeEvolver(alg, dt, nstep, gate, tol, verbosity, state)
 end
 
 #=
@@ -190,41 +190,43 @@ function su_iter(
 end
 
 function Base.iterate(it::TimeEvolver{<:SimpleUpdate}, state = it.state)
-    alg = it.alg
-    iter, t, diff = state.iter, state.t, state.diff
-    check_convergence = (it.tol > 0)
-    if check_convergence
-        if diff < it.tol
-            @info "SU: bond weights have converged."
-            return nothing
-        elseif iter == it.nstep
-            @warn "SU: bond weights have not converged."
-            return nothing
+    return LoggingExtras.withlevel(; it.verbosity) do
+        alg = it.alg
+        iter, t, diff = state.iter, state.t, state.diff
+        check_convergence = (it.tol > 0)
+        if check_convergence
+            if diff < it.tol
+                @infov 1 "SU: bond weights have converged."
+                return nothing
+            elseif iter == it.nstep
+                @warn "SU: bond weights have not converged."
+                return nothing
+            end
+        else
+            (iter == it.nstep) && return nothing
         end
-    else
-        (iter == it.nstep) && return nothing
-    end
-    # perform one iteration and record time
-    time0 = time()
-    ψ, env, ϵ = su_iter(state.ψ, it.gate, it.alg, state.env)
-    diff = compare_weights(env, state.env)
-    elapsed_time = time() - time0
-    # update internal state
-    iter += 1
-    t += it.dt
-    it.state = SUState(iter, t, diff, ψ, env)
-    # output information
-    converged = check_convergence ? (diff < it.tol) : false
-    showinfo = (iter % alg.check_interval == 0) || (iter == 1) || (iter == it.nstep) || converged
-    if showinfo
-        @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))"
-        @info @sprintf(
-            "SU iter %-7d: dt = %.0e, |Δλ| = %.3e. Time = %.3f s/it",
-            iter, it.dt, diff, elapsed_time
-        )
+        # perform one iteration and record time
+        time0 = time()
+        ψ, env, ϵ = su_iter(state.ψ, it.gate, it.alg, state.env)
+        diff = compare_weights(env, state.env)
+        elapsed_time = time() - time0
+        # update internal state
+        iter += 1
+        t += it.dt
+        it.state = SUState(iter, t, diff, ψ, env)
+        # output information
+        converged = check_convergence ? (diff < it.tol) : false
+        showinfo = (iter % alg.check_interval == 0) || (iter == 1) || (iter == it.nstep) || converged
+        if showinfo
+            @infov 1 "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))"
+            @infov 1 @sprintf(
+                "SU iter %-7d: dt = %.0e, |Δλ| = %.3e. Time = %.3f s/it",
+                iter, it.dt, diff, elapsed_time
+            )
+        end
+        info = (; t, ϵ)
+        return (ψ, env, info), it.state
     end
-    info = (; t, ϵ)
-    return (ψ, env, info), it.state
 end
 
 """
@@ -276,7 +278,7 @@ end
     time_evolve(
         ψ₀::Union{InfinitePEPS, InfinitePEPO}, H::LocalOperator, 
         dt::Float64, nstep::Int, alg::SimpleUpdate, env₀::SUWeight;
-        tol::Float64 = 0.0, t₀::Float64 = 0.0
+        tol::Float64 = 0.0, t₀::Float64 = 0.0, verbosity::Int = 1
     ) -> (ψ, env, info)
 
 Perform time evolution on the initial state `ψ₀` and initial environment `env₀`
@@ -292,8 +294,8 @@ with Hamiltonian `H`, using SimpleUpdate algorithm `alg`, time step `dt` for
 function MPSKit.time_evolve(
         ψ₀::InfiniteState, H::LocalOperator, dt::Float64, nstep::Int,
         alg::SimpleUpdate, env₀::SUWeight;
-        tol::Float64 = 0.0, t₀::Float64 = 0.0
+        tol::Float64 = 0.0, t₀::Float64 = 0.0, verbosity::Int = 1
     )
-    it = TimeEvolver(ψ₀, H, dt, nstep, alg, env₀; tol, t₀)
+    it = TimeEvolver(ψ₀, H, dt, nstep, alg, env₀; tol, t₀, verbosity)
     return time_evolve(it)
 end
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index cb4ae3ce..5e75bdf2 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -25,6 +25,8 @@ mutable struct TimeEvolver{TE <: TimeEvolution, G, S}
     gate::G
     # Convergence tolerance (change of weight or energy from last iteration)
     tol::Float64
+    # verbosity level for showing information during evolution
+    verbosity::Int
     # PEPS/PEPO (and environment)
     state::S
 end
diff --git a/test/timeevol/sitedep_truncation.jl b/test/timeevol/sitedep_truncation.jl
index fbc86fb8..3207ac64 100644
--- a/test/timeevol/sitedep_truncation.jl
+++ b/test/timeevol/sitedep_truncation.jl
@@ -28,7 +28,7 @@ end
     bonddims = stack([[6 4; 4 6], [5 7; 7 5]]; dims = 1)
     trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims))
     alg = SimpleUpdate(; trunc, bipartite = true)
-    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0; verbosity = 0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
     # check bipartite structure is preserved
@@ -55,12 +55,12 @@ end
     trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims))
     # 2-site SU
     alg = SimpleUpdate(; trunc)
-    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0; verbosity = 0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
     # 3-site SU
     alg = SimpleUpdate(; trunc, force_3site = true)
-    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0; verbosity = 0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
 end

From b125ed090d6fe373985d51d37f9e6ab7d5a5f1e7 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Mon, 10 Nov 2025 09:33:04 +0800
Subject: [PATCH 19/36] Store Hamiltonian in TimeEvolver

---
 src/algorithms/time_evolution/simpleupdate.jl |  9 ++++-----
 src/algorithms/time_evolution/time_evolve.jl  | 12 +++++++-----
 2 files changed, 11 insertions(+), 10 deletions(-)

diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index cdf73fb6..219718b5 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -29,7 +29,7 @@ struct SUState{S}
     iter::Int
     # evolved time
     t::Float64
-    # measure of difference (of SUWeight, energy, etc.) from last iteration
+    # SUWeight difference from last iteration
     diff::Float64
     # PEPS/PEPO
     ψ::S
@@ -57,7 +57,7 @@ function TimeEvolver(
         tol::Float64 = 0.0, t₀::Float64 = 0.0, verbosity::Int = 1
     )
     # sanity checks
-    _timeevol_sanity_check(ψ₀, physicalspace(H), tol, alg)
+    _timeevol_sanity_check(ψ₀, H, tol, alg)
     # create Trotter gates
     nnonly = is_nearest_neighbour(H)
     use_3site = alg.force_3site || !nnonly
@@ -76,7 +76,7 @@ function TimeEvolver(
         get_expham(H, dt′)
     end
     state = SUState(0, t₀, NaN, ψ₀, env₀)
-    return TimeEvolver(alg, dt, nstep, gate, tol, verbosity, state)
+    return TimeEvolver(alg, dt, nstep, H, gate, tol, verbosity, state)
 end
 
 #=
@@ -245,8 +245,7 @@ function MPSKit.timestep(
         it::TimeEvolver{<:SimpleUpdate}, ψ::InfiniteState, env::SUWeight;
         iter::Int = it.state.iter, t::Float64 = it.state.t
     )
-    Pspaces = physicalspace(it.state.ψ)
-    _timeevol_sanity_check(ψ, Pspaces, it.tol, it.alg)
+    _timeevol_sanity_check(ψ, it.ham, it.tol, it.alg)
     state = SUState(iter, t, NaN, ψ, env)
     result = iterate(it, state)
     if result === nothing
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index 5e75bdf2..b53b3bd5 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -6,7 +6,7 @@ Abstract super type for time evolution algorithms of InfinitePEPS or InfinitePEP
 abstract type TimeEvolution end
 
 """
-    mutable struct TimeEvolver{TE <: TimeEvolution, G, S}
+    mutable struct TimeEvolver{TE <: TimeEvolution, H, G, S}
 
 Iterator for Trotter-based time evolution of InfinitePEPS or InfinitePEPO.
 
@@ -14,13 +14,15 @@ Iterator for Trotter-based time evolution of InfinitePEPS or InfinitePEPO.
 
 $(TYPEDFIELDS)
 """
-mutable struct TimeEvolver{TE <: TimeEvolution, G, S}
+mutable struct TimeEvolver{TE <: TimeEvolution, H, G, S}
     # Time evolution algorithm
     alg::TE
     # Trotter time step
     dt::Float64
     # Maximal iteration steps
     nstep::Int
+    # Hamiltonian
+    ham::H
     # Trotter gates
     gate::G
     # Convergence tolerance (change of weight or energy from last iteration)
@@ -34,11 +36,11 @@ end
 Base.iterate(it::TimeEvolver) = iterate(it, it.state)
 
 function _timeevol_sanity_check(
-        ψ₀::InfiniteState, Pspaces::M, tol::Float64, alg::A
-    ) where {A <: TimeEvolution, M <: AbstractMatrix{<:ElementarySpace}}
+        ψ₀::InfiniteState, ham::LocalOperator, tol::Float64, alg::A
+    ) where {A <: TimeEvolution}
     Nr, Nc, = size(ψ₀)
     @assert (Nr >= 2 && Nc >= 2) "Unit cell size for simple update should be no smaller than (2, 2)."
-    @assert Pspaces == physicalspace(ψ₀) "Physical spaces of `ψ₀` do not match `Pspaces`."
+    @assert physicalspace(ham) == physicalspace(ψ₀) "Physical spaces of `ψ₀` do not match `Pspaces`."
     @assert tol >= 0
     if tol > 0
         @assert alg.imaginary_time "`tol` should be 0 for real time evolution."

From c6baaed06a8afe5f3c79fe3afde78014e4e839ae Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 11 Nov 2025 20:22:40 +0800
Subject: [PATCH 20/36] Move convergence check out of TimeEvolver

---
 examples/heisenberg_su/main.jl                |   4 +-
 examples/hubbard_su/main.jl                   |   4 +-
 examples/j1j2_su/main.jl                      |   4 +-
 src/algorithms/time_evolution/simpleupdate.jl | 142 +++++++++---------
 src/algorithms/time_evolution/time_evolve.jl  |  15 +-
 test/examples/j1j2_finiteT.jl                 |   8 +-
 test/examples/tf_ising_finiteT.jl             |   2 +-
 7 files changed, 85 insertions(+), 94 deletions(-)

diff --git a/examples/heisenberg_su/main.jl b/examples/heisenberg_su/main.jl
index 1fcee9e4..d028c459 100644
--- a/examples/heisenberg_su/main.jl
+++ b/examples/heisenberg_su/main.jl
@@ -72,9 +72,9 @@ dts = [1.0e-2, 1.0e-3, 4.0e-4]
 tols = [1.0e-6, 1.0e-8, 1.0e-8]
 nstep = 10000
 trunc_peps = truncerror(; atol = 1.0e-10) & truncrank(Dbond)
-alg = SimpleUpdate(; trunc = trunc_peps, bipartite = true, check_interval = 500)
+alg = SimpleUpdate(; trunc = trunc_peps, bipartite = true)
 for (dt, tol) in zip(dts, tols)
-    global peps, wts, = time_evolve(peps, H, dt, nstep, alg, wts; tol)
+    global peps, wts, = time_evolve(peps, H, dt, nstep, alg, wts; tol, check_interval = 500)
 end
 
 md"""
diff --git a/examples/hubbard_su/main.jl b/examples/hubbard_su/main.jl
index bfb5350a..77f13338 100644
--- a/examples/hubbard_su/main.jl
+++ b/examples/hubbard_su/main.jl
@@ -66,8 +66,8 @@ maxiter = 20000
 
 for (dt, tol, Dbond) in zip(dts, tols, Ds)
     trunc = truncerror(; atol = 1.0e-10) & truncrank(Dbond)
-    alg = SimpleUpdate(; trunc, bipartite = false, check_interval = 2000)
-    global peps, wts, = time_evolve(peps, H, dt, maxiter, alg, wts; tol)
+    alg = SimpleUpdate(; trunc, bipartite = false)
+    global peps, wts, = time_evolve(peps, H, dt, maxiter, alg, wts; tol, check_interval = 2000)
 end
 
 md"""
diff --git a/examples/j1j2_su/main.jl b/examples/j1j2_su/main.jl
index 6b3ec575..e027f5bc 100644
--- a/examples/j1j2_su/main.jl
+++ b/examples/j1j2_su/main.jl
@@ -51,13 +51,13 @@ on the previously evolved PEPS:
 dt, tol, nstep = 1.0e-2, 1.0e-8, 30000
 check_interval = 4000
 trunc_peps = truncerror(; atol = 1.0e-10) & truncrank(Dbond)
-alg = SimpleUpdate(; trunc = trunc_peps, check_interval)
+alg = SimpleUpdate(; trunc = trunc_peps)
 for J2 in 0.1:0.1:0.5
     ## convert Hamiltonian `LocalOperator` to real floats
     H = real(
         j1_j2_model(ComplexF64, symm, InfiniteSquare(Nr, Nc); J1, J2, sublattice = false),
     )
-    global peps, wts, = time_evolve(peps, H, dt, nstep, alg, wts; tol)
+    global peps, wts, = time_evolve(peps, H, dt, nstep, alg, wts; tol, check_interval)
 end
 
 md"""
diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index 219718b5..5befa8bf 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -12,8 +12,6 @@ $(TYPEDFIELDS)
     trunc::TruncationStrategy
     # Switch for imaginary or real time
     imaginary_time::Bool = true
-    # controls interval to print information
-    check_interval::Int = 500
     # force usage of 3-site simple update
     force_3site::Bool = false
     # ---- only applicable to ground state search ----
@@ -29,8 +27,6 @@ struct SUState{S}
     iter::Int
     # evolved time
     t::Float64
-    # SUWeight difference from last iteration
-    diff::Float64
     # PEPS/PEPO
     ψ::S
     # SUWeight environment
@@ -40,24 +36,20 @@ end
 """
     TimeEvolver(
         ψ₀::InfiniteState, H::LocalOperator, dt::Float64, nstep::Int, 
-        alg::SimpleUpdate, env₀::SUWeight; 
-        tol::Float64 = 0.0, t₀::Float64 = 0.0, verbosity = 1
+        alg::SimpleUpdate, env₀::SUWeight; t₀::Float64 = 0.0
     )
 
 Initialize a TimeEvolver with Hamiltonian `H` and simple update `alg`, 
 starting from the initial state `ψ₀` and SUWeight environment `env₀`.
 
 - The initial (real or imaginary) time is specified by `t₀`.
-- Setting `tol > 0` enables convergence check (for imaginary time evolution of InfinitePEPS only).
-    For other usages it should not be changed.
 """
 function TimeEvolver(
         ψ₀::InfiniteState, H::LocalOperator, dt::Float64, nstep::Int,
-        alg::SimpleUpdate, env₀::SUWeight;
-        tol::Float64 = 0.0, t₀::Float64 = 0.0, verbosity::Int = 1
+        alg::SimpleUpdate, env₀::SUWeight; t₀::Float64 = 0.0
     )
     # sanity checks
-    _timeevol_sanity_check(ψ₀, H, tol, alg)
+    _timeevol_sanity_check(ψ₀, H, alg)
     # create Trotter gates
     nnonly = is_nearest_neighbour(H)
     use_3site = alg.force_3site || !nnonly
@@ -75,11 +67,11 @@ function TimeEvolver(
     else
         get_expham(H, dt′)
     end
-    state = SUState(0, t₀, NaN, ψ₀, env₀)
-    return TimeEvolver(alg, dt, nstep, H, gate, tol, verbosity, state)
+    state = SUState(0, t₀, ψ₀, env₀)
+    return TimeEvolver(alg, dt, nstep, gate, state)
 end
 
-#=
+"""
 Simple update of the x-bond between `[r,c]` and `[r,c+1]`.
 
 ```
@@ -87,7 +79,7 @@ Simple update of the x-bond between `[r,c]` and `[r,c+1]`.
     -- T[r,c] -- T[r,c+1] --
         |           |
 ```
-=#
+"""
 function _su_xbond!(
         state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
         row::Int, col::Int, trunc::TruncationStrategy; gate_ax::Int = 1
@@ -116,7 +108,7 @@ function _su_xbond!(
     return ϵ
 end
 
-#=
+"""
 Simple update of the y-bond between `[r,c]` and `[r-1,c]`.
 ```
         |
@@ -125,7 +117,7 @@ Simple update of the y-bond between `[r,c]` and `[r-1,c]`.
     -- T[r,c] ---
         |
 ```
-=#
+"""
 function _su_ybond!(
         state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
         row::Int, col::Int, trunc::TruncationStrategy; gate_ax::Int = 1
@@ -154,6 +146,9 @@ function _su_ybond!(
     return ϵ
 end
 
+"""
+One iteration of simple update
+"""
 function su_iter(
         state::InfiniteState, gate::LocalOperator, alg::SimpleUpdate, env::SUWeight
     )
@@ -190,43 +185,15 @@ function su_iter(
 end
 
 function Base.iterate(it::TimeEvolver{<:SimpleUpdate}, state = it.state)
-    return LoggingExtras.withlevel(; it.verbosity) do
-        alg = it.alg
-        iter, t, diff = state.iter, state.t, state.diff
-        check_convergence = (it.tol > 0)
-        if check_convergence
-            if diff < it.tol
-                @infov 1 "SU: bond weights have converged."
-                return nothing
-            elseif iter == it.nstep
-                @warn "SU: bond weights have not converged."
-                return nothing
-            end
-        else
-            (iter == it.nstep) && return nothing
-        end
-        # perform one iteration and record time
-        time0 = time()
-        ψ, env, ϵ = su_iter(state.ψ, it.gate, it.alg, state.env)
-        diff = compare_weights(env, state.env)
-        elapsed_time = time() - time0
-        # update internal state
-        iter += 1
-        t += it.dt
-        it.state = SUState(iter, t, diff, ψ, env)
-        # output information
-        converged = check_convergence ? (diff < it.tol) : false
-        showinfo = (iter % alg.check_interval == 0) || (iter == 1) || (iter == it.nstep) || converged
-        if showinfo
-            @infov 1 "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))"
-            @infov 1 @sprintf(
-                "SU iter %-7d: dt = %.0e, |Δλ| = %.3e. Time = %.3f s/it",
-                iter, it.dt, diff, elapsed_time
-            )
-        end
-        info = (; t, ϵ)
-        return (ψ, env, info), it.state
-    end
+    iter, t = state.iter, state.t
+    (iter == it.nstep) && return nothing
+    ψ, env, ϵ = su_iter(state.ψ, it.gate, it.alg, state.env)
+    # update internal state
+    iter += 1
+    t += it.dt
+    it.state = SUState(iter, t, ψ, env)
+    info = (; t, ϵ)
+    return (ψ, env, info), it.state
 end
 
 """
@@ -245,8 +212,8 @@ function MPSKit.timestep(
         it::TimeEvolver{<:SimpleUpdate}, ψ::InfiniteState, env::SUWeight;
         iter::Int = it.state.iter, t::Float64 = it.state.t
     )
-    _timeevol_sanity_check(ψ, it.ham, it.tol, it.alg)
-    state = SUState(iter, t, NaN, ψ, env)
+    _timeevol_sanity_check(ψ, it.ham, it.alg)
+    state = SUState(iter, t, ψ, env)
     result = iterate(it, state)
     if result === nothing
         @warn "TimeEvolver `it` has already reached the end."
@@ -257,27 +224,62 @@ function MPSKit.timestep(
 end
 
 """
-    time_evolve(it::TimeEvolver{<:SimpleUpdate})
+    time_evolve(
+        it::TimeEvolver{<:SimpleUpdate}; 
+        tol::Float64 = 0.0, check_interval::Int = 500
+    )
 
-Perform time evolution until convergence or the set number of iterations
-using the specified TimeEvolver iterator `it` directly.
+Perform time evolution to the end of TimeEvolver iterator `it`,
+or until convergence of SUWeight set by a positive `tol`.
 """
-function MPSKit.time_evolve(it::TimeEvolver{<:SimpleUpdate})
+function MPSKit.time_evolve(
+        it::TimeEvolver{<:SimpleUpdate};
+        tol::Float64 = 0.0, check_interval::Int = 500
+    )
     time_start = time()
-    result = nothing
-    for state in it
-        result = state
+    check_convergence = (tol > 0)
+    if check_convergence
+        @assert (it.state.ψ isa InfinitePEPS) && it.alg.imaginary_time
+    end
+    env0, time0 = it.state.env, time()
+    for (ψ, env, info) in it
+        iter = it.state.iter
+        diff = compare_weights(env0, env)
+        stop = (iter == it.nstep) || (diff < tol)
+        showinfo = (iter % check_interval == 0) || (iter == 1) || stop
+        time1 = time()
+        if showinfo
+            @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))"
+            @info @sprintf(
+                "SU iter %-7d: dt = %.0e, |Δλ| = %.3e. Time = %.3f s/it",
+                iter, it.dt, diff, time1 - time0
+            )
+        end
+        if check_convergence
+            if (iter == it.nstep) && (diff >= tol)
+                @warn "SU: bond weights have not converged."
+            end
+            if diff < tol
+                @info "SU: bond weights have converged."
+            end
+        end
+        if stop
+            time_end = time()
+            @info @sprintf("Simple update finished. Total time elasped: %.2f s", time_end - time_start)
+            return ψ, env, info
+        else
+            env0 = env
+        end
+        time0 = time()
     end
-    time_end = time()
-    @info @sprintf("Simple update finished. Total time elasped: %.2f s", time_end - time_start)
-    return result
+    return
 end
 
 """
     time_evolve(
         ψ₀::Union{InfinitePEPS, InfinitePEPO}, H::LocalOperator, 
         dt::Float64, nstep::Int, alg::SimpleUpdate, env₀::SUWeight;
-        tol::Float64 = 0.0, t₀::Float64 = 0.0, verbosity::Int = 1
+        tol::Float64 = 0.0, t₀::Float64 = 0.0, check_interval::Int = 500
     ) -> (ψ, env, info)
 
 Perform time evolution on the initial state `ψ₀` and initial environment `env₀`
@@ -293,8 +295,8 @@ with Hamiltonian `H`, using SimpleUpdate algorithm `alg`, time step `dt` for
 function MPSKit.time_evolve(
         ψ₀::InfiniteState, H::LocalOperator, dt::Float64, nstep::Int,
         alg::SimpleUpdate, env₀::SUWeight;
-        tol::Float64 = 0.0, t₀::Float64 = 0.0, verbosity::Int = 1
+        tol::Float64 = 0.0, t₀::Float64 = 0.0, check_interval::Int = 500
     )
-    it = TimeEvolver(ψ₀, H, dt, nstep, alg, env₀; tol, t₀, verbosity)
-    return time_evolve(it)
+    it = TimeEvolver(ψ₀, H, dt, nstep, alg, env₀; t₀)
+    return time_evolve(it; tol, check_interval)
 end
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index b53b3bd5..a14582fb 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -14,21 +14,15 @@ Iterator for Trotter-based time evolution of InfinitePEPS or InfinitePEPO.
 
 $(TYPEDFIELDS)
 """
-mutable struct TimeEvolver{TE <: TimeEvolution, H, G, S}
+mutable struct TimeEvolver{TE <: TimeEvolution, G, S}
     # Time evolution algorithm
     alg::TE
     # Trotter time step
     dt::Float64
     # Maximal iteration steps
     nstep::Int
-    # Hamiltonian
-    ham::H
     # Trotter gates
     gate::G
-    # Convergence tolerance (change of weight or energy from last iteration)
-    tol::Float64
-    # verbosity level for showing information during evolution
-    verbosity::Int
     # PEPS/PEPO (and environment)
     state::S
 end
@@ -36,16 +30,11 @@ end
 Base.iterate(it::TimeEvolver) = iterate(it, it.state)
 
 function _timeevol_sanity_check(
-        ψ₀::InfiniteState, ham::LocalOperator, tol::Float64, alg::A
+        ψ₀::InfiniteState, ham::LocalOperator, alg::A
     ) where {A <: TimeEvolution}
     Nr, Nc, = size(ψ₀)
     @assert (Nr >= 2 && Nc >= 2) "Unit cell size for simple update should be no smaller than (2, 2)."
     @assert physicalspace(ham) == physicalspace(ψ₀) "Physical spaces of `ψ₀` do not match `Pspaces`."
-    @assert tol >= 0
-    if tol > 0
-        @assert alg.imaginary_time "`tol` should be 0 for real time evolution."
-        @assert ψ₀ isa InfinitePEPS "`tol` should be 0 for time evolution of InfinitePEPO."
-    end
     if hasfield(typeof(alg), :gate_bothsides) && alg.gate_bothsides
         @assert ψ₀ isa InfinitePEPO "alg.gate_bothsides = true is only compatible with PEPO."
     end
diff --git a/test/examples/j1j2_finiteT.jl b/test/examples/j1j2_finiteT.jl
index 8ffa3836..1ec19bfb 100644
--- a/test/examples/j1j2_finiteT.jl
+++ b/test/examples/j1j2_finiteT.jl
@@ -32,9 +32,9 @@ check_interval = 100
 
 # PEPO approach
 dt, nstep = 1.0e-3, 600
-alg = SimpleUpdate(; trunc = trunc_pepo, gate_bothsides = true, check_interval)
+alg = SimpleUpdate(; trunc = trunc_pepo, gate_bothsides = true)
 evolver = TimeEvolver(pepo0, ham, dt, nstep, alg, wts0)
-pepo, wts, info = time_evolve(evolver)
+pepo, wts, info = time_evolve(evolver; check_interval)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 energy = expectation_value(pepo, ham, env) / (Nr * Nc)
 @info "β = $(dt * nstep): tr(ρH) = $(energy)"
@@ -42,9 +42,9 @@ energy = expectation_value(pepo, ham, env) / (Nr * Nc)
 @test energy ≈ bm[2] atol = 5.0e-3
 
 # PEPS (purified PEPO) approach
-alg = SimpleUpdate(; trunc = trunc_pepo, gate_bothsides = false, check_interval)
+alg = SimpleUpdate(; trunc = trunc_pepo, gate_bothsides = false)
 evolver = TimeEvolver(pepo0, ham, dt, nstep, alg, wts0)
-pepo, wts, = time_evolve(evolver)
+pepo, wts, = time_evolve(evolver; check_interval)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 energy = expectation_value(pepo, ham, env) / (Nr * Nc)
 @info "β = $(dt * nstep) / 2: tr(ρH) = $(energy)"
diff --git a/test/examples/tf_ising_finiteT.jl b/test/examples/tf_ising_finiteT.jl
index 766fb632..109b4109 100644
--- a/test/examples/tf_ising_finiteT.jl
+++ b/test/examples/tf_ising_finiteT.jl
@@ -25,7 +25,7 @@ end
 
 function converge_env(state, χ::Int)
     trunc1 = truncrank(4) & truncerror(; atol = 1.0e-12)
-    env0 = CTMRGEnv(randn, Float64, state, ℂ^4)
+    env0 = CTMRGEnv(ones, Float64, state, ℂ^1)
     env, = leading_boundary(env0, state; alg = :sequential, trunc = trunc1, tol = 1.0e-10)
     trunc2 = truncrank(χ) & truncerror(; atol = 1.0e-12)
     env, = leading_boundary(env, state; alg = :sequential, trunc = trunc2, tol = 1.0e-10)

From e841bf147d4f54a07ad016a5a314a427c752efa9 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 11 Nov 2025 20:27:06 +0800
Subject: [PATCH 21/36] Restore docstrings

---
 .../time_evolution/simpleupdate3site.jl       | 32 +++++++++----------
 1 file changed, 16 insertions(+), 16 deletions(-)

diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl
index c75663b9..9be45a7f 100644
--- a/src/algorithms/time_evolution/simpleupdate3site.jl
+++ b/src/algorithms/time_evolution/simpleupdate3site.jl
@@ -124,14 +124,14 @@ Then the fidelity is just
     F(ψ̃) = (norm(s̃[n], 2) / norm(s[n], 2))^2
 ```
 =#
-#=
+"""
 Perform QR decomposition through a PEPS tensor
 ```
              ╱           ╱
     --R0----M---  →  ---Q--*-R1--
           ╱ |         ╱ |
 ```
-=#
+"""
 function qr_through(
         R0::MPSBondTensor, M::GenericMPSTensor{S, 4}; normalize::Bool = true
     ) where {S <: ElementarySpace}
@@ -154,14 +154,14 @@ function qr_through(
     return r
 end
 
-#=
+"""
 Perform LQ decomposition through a tensor
 ```
              ╱           ╱
     --L0-*--Q---  ←  ---M--*-L1--
           ╱ |         ╱ |
 ```
-=#
+"""
 function lq_through(
         M::GenericMPSTensor{S, 4}, L1::MPSBondTensor; normalize::Bool = true
     ) where {S <: ElementarySpace}
@@ -186,9 +186,9 @@ function lq_through(
     return l
 end
 
-#=
+"""
 Given a cluster `Ms`, find all `R`, `L` matrices on each internal bond
-=#
+"""
 function _get_allRLs(Ms::Vector{T}) where {T <: GenericMPSTensor{<:ElementarySpace, 4}}
     # M1 -- (R1,L1) -- M2 -- (R2,L2) -- M3
     N = length(Ms)
@@ -207,7 +207,7 @@ function _get_allRLs(Ms::Vector{T}) where {T <: GenericMPSTensor{<:ElementarySpa
     return Rs, Ls
 end
 
-#=
+"""
 Given the tensors `R`, `L` on a bond, construct 
 the projectors `Pa`, `Pb` and the new bond weight `s`
 such that the contraction of `Pa`, `s`, `Pb` is identity when `trunc = notrunc`,
@@ -220,7 +220,7 @@ The arrows between `Pa`, `s`, `Pb` are
     rev = true:  - Pa --→-- Pb - 
                     2 → s → 1
 ```
-=#
+"""
 function _proj_from_RL(
         r::MPSBondTensor, l::MPSBondTensor;
         trunc::TruncationStrategy = notrunc(), rev::Bool = false,
@@ -240,10 +240,10 @@ function _proj_from_RL(
     return Pa, s, Pb, ϵ
 end
 
-#=
+"""
 Given a cluster `Ms` and the pre-calculated `R`, `L` bond matrices,
 find all projectors `Pa`, `Pb` and Schmidt weights `wts` on internal bonds.
-=#
+"""
 function _get_allprojs(
         Ms, Rs, Ls, truncs::Vector{E}, revs::Vector{Bool}
     ) where {E <: TruncationStrategy}
@@ -266,9 +266,9 @@ function _get_allprojs(
     return Pas, Pbs, wts, ϵs
 end
 
-#=
+"""
 Find projectors to truncate internal bonds of the cluster `Ms`.
-=#
+"""
 function _cluster_truncate!(
         Ms::Vector{T}, truncs::Vector{E}, revs::Vector{Bool}
     ) where {T <: GenericMPSTensor{<:ElementarySpace, 4}, E <: TruncationStrategy}
@@ -283,7 +283,7 @@ function _cluster_truncate!(
     return wts, ϵs, Pas, Pbs
 end
 
-#=
+"""
 Apply the gate MPO `gs` on the cluster `Ms`.
 When `gate_ax` is 1 or 2, the gate acts from the physical codomain or domain side.
 
@@ -307,7 +307,7 @@ In the cluster, the axes of each tensor use the MPS order
       4  2            5  2
     M[1 2 3 4; 5]  M[1 2 3 4 5; 6]
 ```
-=#
+"""
 function _apply_gatempo!(
         Ms::Vector{T1}, gs::Vector{T2}; gate_ax::Int = 1
     ) where {T1 <: GenericMPSTensor{<:ElementarySpace, 4}, T2 <: AbstractTensorMap}
@@ -437,7 +437,7 @@ const perms_se_pepo = map(invperms_se_pepo) do (p1, p2)
     p = invperm((p1..., p2...))
     return (p[1:(end - 1)], (p[end],))
 end
-#=
+"""
 Obtain the 3-site cluster in the "southeast corner" of a square plaquette.
 ``` 
     r-1         M3
@@ -446,7 +446,7 @@ Obtain the 3-site cluster in the "southeast corner" of a square plaquette.
     r   M1 -←- M2
         c      c+1
 ```
-=#
+"""
 function get_3site_se(state::InfiniteState, env::SUWeight, row::Int, col::Int)
     Nr, Nc = size(state)
     rm1, cp1 = _prev(row, Nr), _next(col, Nc)

From ef5f353f9d1b031fc363eceb528241b116648fa5 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 11 Nov 2025 20:43:28 +0800
Subject: [PATCH 22/36] Fix `timestep`

---
 src/algorithms/time_evolution/simpleupdate.jl | 4 ++--
 src/algorithms/time_evolution/time_evolve.jl  | 6 +++---
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index 5befa8bf..4046b8ee 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -49,7 +49,7 @@ function TimeEvolver(
         alg::SimpleUpdate, env₀::SUWeight; t₀::Float64 = 0.0
     )
     # sanity checks
-    _timeevol_sanity_check(ψ₀, H, alg)
+    _timeevol_sanity_check(ψ₀, physicalspace(H), alg)
     # create Trotter gates
     nnonly = is_nearest_neighbour(H)
     use_3site = alg.force_3site || !nnonly
@@ -212,7 +212,7 @@ function MPSKit.timestep(
         it::TimeEvolver{<:SimpleUpdate}, ψ::InfiniteState, env::SUWeight;
         iter::Int = it.state.iter, t::Float64 = it.state.t
     )
-    _timeevol_sanity_check(ψ, it.ham, it.alg)
+    _timeevol_sanity_check(ψ, physicalspace(it.state.ψ), it.alg)
     state = SUState(iter, t, ψ, env)
     result = iterate(it, state)
     if result === nothing
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index a14582fb..d1fb7eb0 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -30,11 +30,11 @@ end
 Base.iterate(it::TimeEvolver) = iterate(it, it.state)
 
 function _timeevol_sanity_check(
-        ψ₀::InfiniteState, ham::LocalOperator, alg::A
-    ) where {A <: TimeEvolution}
+        ψ₀::InfiniteState, Pspaces::M, alg::A
+    ) where {A <: TimeEvolution, M <: AbstractMatrix{<:ElementarySpace}}
     Nr, Nc, = size(ψ₀)
     @assert (Nr >= 2 && Nc >= 2) "Unit cell size for simple update should be no smaller than (2, 2)."
-    @assert physicalspace(ham) == physicalspace(ψ₀) "Physical spaces of `ψ₀` do not match `Pspaces`."
+    @assert Pspaces == physicalspace(ψ₀) "Physical spaces of `ψ₀` do not match `Pspaces`."
     if hasfield(typeof(alg), :gate_bothsides) && alg.gate_bothsides
         @assert ψ₀ isa InfinitePEPO "alg.gate_bothsides = true is only compatible with PEPO."
     end

From 5ec704ec1e417a40be5288357fd503c42580a49b Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 11 Nov 2025 20:58:27 +0800
Subject: [PATCH 23/36] Fix tests

---
 test/bondenv/benv_fu.jl             | 6 ++----
 test/timeevol/cluster_projectors.jl | 8 ++++----
 2 files changed, 6 insertions(+), 8 deletions(-)

diff --git a/test/bondenv/benv_fu.jl b/test/bondenv/benv_fu.jl
index 8d74f96a..9212f3b1 100644
--- a/test/bondenv/benv_fu.jl
+++ b/test/bondenv/benv_fu.jl
@@ -13,10 +13,8 @@ function get_hubbard_state(t::Float64 = 1.0, U::Float64 = 8.0)
     Vphy = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 2, (1, 1 // 2) => 1, (1, -1 // 2) => 1)
     peps = InfinitePEPS(rand, ComplexF64, Vphy, Vphy; unitcell = (Nr, Nc))
     wts = SUWeight(peps)
-    alg = SimpleUpdate(;
-        trunc = truncerror(; atol = 1.0e-10) & truncrank(4), check_interval = 2000
-    )
-    evolver = TimeEvolver(peps, H, 1.0e-2, 10000, alg, wts; tol = 1.0e-8)
+    alg = SimpleUpdate(; trunc = truncerror(; atol = 1.0e-10) & truncrank(4))
+    evolver = TimeEvolver(peps, H, 1.0e-2, 10000, alg, wts; tol = 1.0e-8, check_interval = 2000)
     peps, = time_evolve(evolver)
     normalize!.(peps.A, Inf)
     return peps
diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl
index 71f4163a..279907c1 100644
--- a/test/timeevol/cluster_projectors.jl
+++ b/test/timeevol/cluster_projectors.jl
@@ -120,8 +120,8 @@ end
     nstep = 5000
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
         trunc = truncerror(; atol = 1.0e-10) & truncrank(n == 1 ? 4 : 2)
-        alg = SimpleUpdate(; trunc, bipartite = true, check_interval = 1000)
-        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts; tol)
+        alg = SimpleUpdate(; trunc, bipartite = true)
+        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts; tol, check_interval = 1000)
     end
     normalize!.(peps.A, Inf)
     env = CTMRGEnv(wts, peps)
@@ -132,9 +132,9 @@ end
     dts = [1.0e-2, 5.0e-3]
     tols = [1.0e-8, 1.0e-8]
     trunc = truncerror(; atol = 1.0e-10) & truncrank(2)
-    alg = SimpleUpdate(; trunc, check_interval = 1000, force_3site = true)
+    alg = SimpleUpdate(; trunc, force_3site = true)
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
-        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts; tol)
+        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts; tol, check_interval = 1000)
     end
     normalize!.(peps.A, Inf)
     env, = leading_boundary(env, peps; tol = ctmrg_tol, trunc = trunc_env)

From 97dd4208c25fdee4dd3404357a01b975ce3de768 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 11 Nov 2025 21:52:35 +0800
Subject: [PATCH 24/36] Fix more tests

---
 test/bondenv/benv_fu.jl             | 4 ++--
 test/timeevol/sitedep_truncation.jl | 6 +++---
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/test/bondenv/benv_fu.jl b/test/bondenv/benv_fu.jl
index 9212f3b1..1a531032 100644
--- a/test/bondenv/benv_fu.jl
+++ b/test/bondenv/benv_fu.jl
@@ -14,8 +14,8 @@ function get_hubbard_state(t::Float64 = 1.0, U::Float64 = 8.0)
     peps = InfinitePEPS(rand, ComplexF64, Vphy, Vphy; unitcell = (Nr, Nc))
     wts = SUWeight(peps)
     alg = SimpleUpdate(; trunc = truncerror(; atol = 1.0e-10) & truncrank(4))
-    evolver = TimeEvolver(peps, H, 1.0e-2, 10000, alg, wts; tol = 1.0e-8, check_interval = 2000)
-    peps, = time_evolve(evolver)
+    evolver = TimeEvolver(peps, H, 1.0e-2, 10000, alg, wts)
+    peps, = time_evolve(evolver; tol = 1.0e-8, check_interval = 2000)
     normalize!.(peps.A, Inf)
     return peps
 end
diff --git a/test/timeevol/sitedep_truncation.jl b/test/timeevol/sitedep_truncation.jl
index 3207ac64..fbc86fb8 100644
--- a/test/timeevol/sitedep_truncation.jl
+++ b/test/timeevol/sitedep_truncation.jl
@@ -28,7 +28,7 @@ end
     bonddims = stack([[6 4; 4 6], [5 7; 7 5]]; dims = 1)
     trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims))
     alg = SimpleUpdate(; trunc, bipartite = true)
-    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0; verbosity = 0)
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
     # check bipartite structure is preserved
@@ -55,12 +55,12 @@ end
     trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims))
     # 2-site SU
     alg = SimpleUpdate(; trunc)
-    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0; verbosity = 0)
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
     # 3-site SU
     alg = SimpleUpdate(; trunc, force_3site = true)
-    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0; verbosity = 0)
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
 end

From 8207b73f16ddfba9edeac5ffc634bd964ba541ed Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Wed, 12 Nov 2025 09:37:01 +0800
Subject: [PATCH 25/36] Finer parameters for SUState

---
 src/algorithms/time_evolution/simpleupdate.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index 4046b8ee..a184d4e8 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -22,7 +22,7 @@ $(TYPEDFIELDS)
 end
 
 # internal state of simple update algorithm
-struct SUState{S}
+struct SUState{S <: InfiniteState, E <: SUWeight}
     # number of performed iterations
     iter::Int
     # evolved time
@@ -30,7 +30,7 @@ struct SUState{S}
     # PEPS/PEPO
     ψ::S
     # SUWeight environment
-    env::SUWeight
+    env::E
 end
 
 """

From 1d3e962a5d89bc3dfa5e7bd380cda6b3ca053fc1 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Wed, 12 Nov 2025 17:32:12 +0800
Subject: [PATCH 26/36] Make timestep test look better

---
 src/algorithms/time_evolution/simpleupdate.jl |  2 +-
 test/timeevol/timestep.jl                     | 34 +++++++++----------
 2 files changed, 18 insertions(+), 18 deletions(-)

diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index a184d4e8..ac922fa5 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -239,7 +239,7 @@ function MPSKit.time_evolve(
     time_start = time()
     check_convergence = (tol > 0)
     if check_convergence
-        @assert (it.state.ψ isa InfinitePEPS) && it.alg.imaginary_time
+        @assert (it.state.ψ isa InfinitePEPS) && it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking."
     end
     env0, time0 = it.state.env, time()
     for (ψ, env, info) in it
diff --git a/test/timeevol/timestep.jl b/test/timeevol/timestep.jl
index ede88db6..c1513bd9 100644
--- a/test/timeevol/timestep.jl
+++ b/test/timeevol/timestep.jl
@@ -3,9 +3,13 @@ using Random
 using TensorKit
 using PEPSKit
 
-function test_timestep(ψ0, env0, H)
-    trunc = truncerror(; atol = 1.0e-10) & truncrank(4)
-    alg = SimpleUpdate(; trunc)
+@testset "SimpleUpdate timestep" begin
+    Nr, Nc = 2, 2
+    H = real(heisenberg_XYZ(ComplexF64, Trivial, InfiniteSquare(Nr, Nc); Jx = 1, Jy = 1, Jz = 1))
+    Pspace, Vspace = ℂ^2, ℂ^4
+    ψ0 = InfinitePEPS(rand, Float64, Pspace, Vspace; unitcell = (Nr, Nc))
+    env0 = SUWeight(ψ0)
+    alg = SimpleUpdate(; trunc = truncerror(; atol = 1.0e-10) & truncrank(4))
     dt, nstep = 1.0e-2, 50
     # manual timestep
     evolver = TimeEvolver(ψ0, H, dt, nstep, alg, env0)
@@ -13,21 +17,17 @@ function test_timestep(ψ0, env0, H)
     for iter in 0:(nstep - 1)
         ψ1, env1, info1 = timestep(evolver, ψ1, env1; iter)
     end
-    @info info1
     # time_evolve
     ψ2, env2, info2 = time_evolve(ψ0, H, dt, nstep, alg, env0)
-    @info info2
+    # for-loop syntax
+    ## manually reset internal state of evolver
+    evolver.state = PEPSKit.SUState(0, 0.0, ψ0, env0)
+    ψ3, env3, info3 = nothing, nothing, nothing
+    for state in evolver
+        ψ3, env3, info3 = state
+    end
     # results should be *exactly* the same
-    @test ψ1 == ψ2
-    @test env1 == env2
-    @test info1 == info2
-    return nothing
+    @test ψ1 == ψ2 == ψ3
+    @test env1 == env2 == env3
+    @test info1 == info2 == info3
 end
-
-Nr, Nc = 2, 2
-H = real(heisenberg_XYZ(ComplexF64, Trivial, InfiniteSquare(Nr, Nc); Jx = 1, Jy = 1, Jz = 1))
-Pspace, Vspace = ℂ^2, ℂ^4
-ψ0 = InfinitePEPS(rand, Float64, Pspace, Vspace; unitcell = (Nr, Nc))
-env0 = SUWeight(ψ0)
-
-test_timestep(ψ0, env0, H)

From c323c798eaaa2faae168e97a6a56abe052559167 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Wed, 12 Nov 2025 23:05:50 +0800
Subject: [PATCH 27/36] Allow complex Trotter time step

---
 src/algorithms/time_evolution/evoltools.jl    |  2 +-
 src/algorithms/time_evolution/simpleupdate.jl | 29 ++++++++++---------
 src/algorithms/time_evolution/time_evolve.jl  |  6 ++--
 3 files changed, 19 insertions(+), 18 deletions(-)

diff --git a/src/algorithms/time_evolution/evoltools.jl b/src/algorithms/time_evolution/evoltools.jl
index d1f66bb5..9659bc95 100644
--- a/src/algorithms/time_evolution/evoltools.jl
+++ b/src/algorithms/time_evolution/evoltools.jl
@@ -1,7 +1,7 @@
 const InfiniteState = Union{InfinitePEPS, InfinitePEPO}
 
 function _get_dt(
-        state::InfiniteState, dt::Float64, imaginary_time::Bool
+        state::InfiniteState, dt::Number, imaginary_time::Bool
     )
     dt′ = (state isa InfinitePEPS) ? dt : (dt / 2)
     if (state isa InfinitePEPO)
diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index ac922fa5..00eafb30 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -10,7 +10,7 @@ $(TYPEDFIELDS)
 @kwdef struct SimpleUpdate <: TimeEvolution
     # Truncation scheme after applying Trotter gates
     trunc::TruncationStrategy
-    # Switch for imaginary or real time
+    # Switch for imaginary time (exp(-H dt)) or real time (exp(-iH dt)) evolution
     imaginary_time::Bool = true
     # force usage of 3-site simple update
     force_3site::Bool = false
@@ -22,11 +22,11 @@ $(TYPEDFIELDS)
 end
 
 # internal state of simple update algorithm
-struct SUState{S <: InfiniteState, E <: SUWeight}
+struct SUState{N <: Number, S <: InfiniteState, E <: SUWeight}
     # number of performed iterations
     iter::Int
     # evolved time
-    t::Float64
+    t::N
     # PEPS/PEPO
     ψ::S
     # SUWeight environment
@@ -35,8 +35,8 @@ end
 
 """
     TimeEvolver(
-        ψ₀::InfiniteState, H::LocalOperator, dt::Float64, nstep::Int, 
-        alg::SimpleUpdate, env₀::SUWeight; t₀::Float64 = 0.0
+        ψ₀::InfiniteState, H::LocalOperator, dt::Number, nstep::Int, 
+        alg::SimpleUpdate, env₀::SUWeight; t₀::Number = 0.0
     )
 
 Initialize a TimeEvolver with Hamiltonian `H` and simple update `alg`, 
@@ -45,8 +45,8 @@ starting from the initial state `ψ₀` and SUWeight environment `env₀`.
 - The initial (real or imaginary) time is specified by `t₀`.
 """
 function TimeEvolver(
-        ψ₀::InfiniteState, H::LocalOperator, dt::Float64, nstep::Int,
-        alg::SimpleUpdate, env₀::SUWeight; t₀::Float64 = 0.0
+        ψ₀::InfiniteState, H::LocalOperator, dt::Number, nstep::Int,
+        alg::SimpleUpdate, env₀::SUWeight; t₀::Number = 0.0
     )
     # sanity checks
     _timeevol_sanity_check(ψ₀, physicalspace(H), alg)
@@ -217,7 +217,7 @@ function MPSKit.timestep(
     result = iterate(it, state)
     if result === nothing
         @warn "TimeEvolver `it` has already reached the end."
-        return ψ, env, (; t, ϵ = NaN)
+        return nothing
     else
         return first(result)
     end
@@ -251,8 +251,9 @@ function MPSKit.time_evolve(
         if showinfo
             @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))"
             @info @sprintf(
-                "SU iter %-7d: dt = %.0e, |Δλ| = %.3e. Time = %.3f s/it",
-                iter, it.dt, diff, time1 - time0
+                "SU iter %-7d: dt = %s, |Δλ| = %.3e. Time = %.3f s/it",
+                # using `string` since `dt` can be complex
+                iter, string(it.dt), diff, time1 - time0
             )
         end
         if check_convergence
@@ -278,8 +279,8 @@ end
 """
     time_evolve(
         ψ₀::Union{InfinitePEPS, InfinitePEPO}, H::LocalOperator, 
-        dt::Float64, nstep::Int, alg::SimpleUpdate, env₀::SUWeight;
-        tol::Float64 = 0.0, t₀::Float64 = 0.0, check_interval::Int = 500
+        dt::Number, nstep::Int, alg::SimpleUpdate, env₀::SUWeight;
+        tol::Float64 = 0.0, t₀::Number = 0.0, check_interval::Int = 500
     ) -> (ψ, env, info)
 
 Perform time evolution on the initial state `ψ₀` and initial environment `env₀`
@@ -293,9 +294,9 @@ with Hamiltonian `H`, using SimpleUpdate algorithm `alg`, time step `dt` for
     including the time evolved since `ψ₀`.
 """
 function MPSKit.time_evolve(
-        ψ₀::InfiniteState, H::LocalOperator, dt::Float64, nstep::Int,
+        ψ₀::InfiniteState, H::LocalOperator, dt::Number, nstep::Int,
         alg::SimpleUpdate, env₀::SUWeight;
-        tol::Float64 = 0.0, t₀::Float64 = 0.0, check_interval::Int = 500
+        tol::Float64 = 0.0, t₀::Number = 0.0, check_interval::Int = 500
     )
     it = TimeEvolver(ψ₀, H, dt, nstep, alg, env₀; t₀)
     return time_evolve(it; tol, check_interval)
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index d1fb7eb0..6b778a20 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -6,7 +6,7 @@ Abstract super type for time evolution algorithms of InfinitePEPS or InfinitePEP
 abstract type TimeEvolution end
 
 """
-    mutable struct TimeEvolver{TE <: TimeEvolution, H, G, S}
+    mutable struct TimeEvolver{TE <: TimeEvolution, N <: Number, G, S}
 
 Iterator for Trotter-based time evolution of InfinitePEPS or InfinitePEPO.
 
@@ -14,11 +14,11 @@ Iterator for Trotter-based time evolution of InfinitePEPS or InfinitePEPO.
 
 $(TYPEDFIELDS)
 """
-mutable struct TimeEvolver{TE <: TimeEvolution, G, S}
+mutable struct TimeEvolver{TE <: TimeEvolution, N <: Number, G, S}
     # Time evolution algorithm
     alg::TE
     # Trotter time step
-    dt::Float64
+    dt::N
     # Maximal iteration steps
     nstep::Int
     # Trotter gates

From f78268d45e61a8bfb800f9753324a0dbcbb23101 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Wed, 12 Nov 2025 23:10:27 +0800
Subject: [PATCH 28/36] Reorder parameters of SUState

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

diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index 00eafb30..aa72b23b 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -22,7 +22,7 @@ $(TYPEDFIELDS)
 end
 
 # internal state of simple update algorithm
-struct SUState{N <: Number, S <: InfiniteState, E <: SUWeight}
+struct SUState{S <: InfiniteState, E <: SUWeight, N <: Number}
     # number of performed iterations
     iter::Int
     # evolved time

From c6dce5c4349a398ee50f9f842b9233b80db9103a Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Thu, 13 Nov 2025 08:47:48 +0800
Subject: [PATCH 29/36] Docstring updates

---
 src/algorithms/time_evolution/simpleupdate.jl | 46 +++++++++++--------
 src/algorithms/time_evolution/time_evolve.jl  | 14 +++---
 src/algorithms/truncation/bond_truncation.jl  |  2 +-
 3 files changed, 36 insertions(+), 26 deletions(-)

diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index aa72b23b..abb54ea1 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -8,28 +8,29 @@ Algorithm struct for simple update (SU) of InfinitePEPS or InfinitePEPO.
 $(TYPEDFIELDS)
 """
 @kwdef struct SimpleUpdate <: TimeEvolution
-    # Truncation scheme after applying Trotter gates
+    "Truncation strategy for bonds updated by Trotter gates"
     trunc::TruncationStrategy
-    # Switch for imaginary time (exp(-H dt)) or real time (exp(-iH dt)) evolution
+    "When true (or false), the Trotter gate is `exp(-H dt)` (or `exp(-iH dt)`)"
     imaginary_time::Bool = true
-    # force usage of 3-site simple update
+    "When true, force the usage of 3-site simple update"
     force_3site::Bool = false
-    # ---- only applicable to ground state search ----
-    # assume bipartite unit cell structure
+    "(Only applicable to InfinitePEPS) When true, assume bipartite unit cell structure"
     bipartite::Bool = false
-    # ---- only applicable to PEPO evolution ----
-    gate_bothsides::Bool = false # when false, purified approach is assumed
+    "(Only applicable to InfinitePEPO) When true (or false), the PEPO is updated as `exp(-H dt/2) * ρ * exp(-H dt/2)` (or `exp(-H dt) ρ`) in each Trotter step"
+    gate_bothsides::Bool = false
 end
 
-# internal state of simple update algorithm
+"""
+Internal state of simple update algorithm
+"""
 struct SUState{S <: InfiniteState, E <: SUWeight, N <: Number}
-    # number of performed iterations
+    "number of performed iterations"
     iter::Int
-    # evolved time
+    "evolved time"
     t::N
-    # PEPS/PEPO
+    "PEPS/PEPO"
     ψ::S
-    # SUWeight environment
+    "SUWeight environment"
     env::E
 end
 
@@ -42,7 +43,7 @@ end
 Initialize a TimeEvolver with Hamiltonian `H` and simple update `alg`, 
 starting from the initial state `ψ₀` and SUWeight environment `env₀`.
 
-- The initial (real or imaginary) time is specified by `t₀`.
+- The initial time is specified by `t₀`.
 """
 function TimeEvolver(
         ψ₀::InfiniteState, H::LocalOperator, dt::Number, nstep::Int,
@@ -73,12 +74,13 @@ end
 
 """
 Simple update of the x-bond between `[r,c]` and `[r,c+1]`.
-
 ```
         |           |
     -- T[r,c] -- T[r,c+1] --
         |           |
 ```
+When `gate_ax = 1` (or `2`), the gate will be applied to 
+the codomain (or domain) physicsl legs of `state`.
 """
 function _su_xbond!(
         state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
@@ -117,6 +119,8 @@ Simple update of the y-bond between `[r,c]` and `[r-1,c]`.
     -- T[r,c] ---
         |
 ```
+When `gate_ax = 1` (or `2`), the gate will be applied to 
+the codomain (or domain) physicsl legs of `state`.
 """
 function _su_ybond!(
         state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
@@ -227,10 +231,14 @@ end
     time_evolve(
         it::TimeEvolver{<:SimpleUpdate}; 
         tol::Float64 = 0.0, check_interval::Int = 500
-    )
+    ) -> (ψ, env, info)
 
 Perform time evolution to the end of TimeEvolver iterator `it`,
 or until convergence of SUWeight set by a positive `tol`.
+
+- Setting `tol > 0` enables convergence check (for imaginary time evolution of InfinitePEPS only).
+    For other usages it should not be changed.
+- `check_interval` sets the number of iterations between outputs of information.
 """
 function MPSKit.time_evolve(
         it::TimeEvolver{<:SimpleUpdate};
@@ -246,7 +254,8 @@ function MPSKit.time_evolve(
         iter = it.state.iter
         diff = compare_weights(env0, env)
         stop = (iter == it.nstep) || (diff < tol)
-        showinfo = (iter % check_interval == 0) || (iter == 1) || stop
+        showinfo = (check_interval > 0) &&
+            ((iter % check_interval == 0) || (iter == 1) || stop)
         time1 = time()
         if showinfo
             @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))"
@@ -289,9 +298,10 @@ with Hamiltonian `H`, using SimpleUpdate algorithm `alg`, time step `dt` for
 
 - Setting `tol > 0` enables convergence check (for imaginary time evolution of InfinitePEPS only).
     For other usages it should not be changed.
-- Using `t₀` to specify the initial (real or imaginary) time of `ψ₀`.
+- Use `t₀` to specify the initial time of `ψ₀`.
+- `check_interval` sets the interval to output information. Output during the evolution can be turned off by setting `check_interval <= 0`.
 - `info` is a NamedTuple containing information of the evolution, 
-    including the time evolved since `ψ₀`.
+    including the time `info.t` evolved since `ψ₀`.
 """
 function MPSKit.time_evolve(
         ψ₀::InfiniteState, H::LocalOperator, dt::Number, nstep::Int,
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index 6b778a20..5071bc1a 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -6,7 +6,7 @@ Abstract super type for time evolution algorithms of InfinitePEPS or InfinitePEP
 abstract type TimeEvolution end
 
 """
-    mutable struct TimeEvolver{TE <: TimeEvolution, N <: Number, G, S}
+    mutable struct TimeEvolver{TE <: TimeEvolution, G, S, N <: Number}
 
 Iterator for Trotter-based time evolution of InfinitePEPS or InfinitePEPO.
 
@@ -14,16 +14,16 @@ Iterator for Trotter-based time evolution of InfinitePEPS or InfinitePEPO.
 
 $(TYPEDFIELDS)
 """
-mutable struct TimeEvolver{TE <: TimeEvolution, N <: Number, G, S}
-    # Time evolution algorithm
+mutable struct TimeEvolver{TE <: TimeEvolution, G, S, N <: Number}
+    "Time evolution algorithm (currently supported: `SimpleUpdate`)"
     alg::TE
-    # Trotter time step
+    "Trotter time step"
     dt::N
-    # Maximal iteration steps
+    "The number of iteration steps"
     nstep::Int
-    # Trotter gates
+    "Trotter gates"
     gate::G
-    # PEPS/PEPO (and environment)
+    "PEPS/PEPO and its environment"
     state::S
 end
 
diff --git a/src/algorithms/truncation/bond_truncation.jl b/src/algorithms/truncation/bond_truncation.jl
index 25fcc190..b2f54ac1 100644
--- a/src/algorithms/truncation/bond_truncation.jl
+++ b/src/algorithms/truncation/bond_truncation.jl
@@ -15,7 +15,7 @@ The truncation algorithm can be constructed from the following keyword arguments
 
 * `trunc::TruncationStrategy`: SVD truncation strategy when initilizing the truncated tensors connected by the bond.
 * `maxiter::Int=50` : Maximal number of ALS iterations.
-* `tol::Float64=1e-15` : ALS converges when fidelity change between two FET iterations is smaller than `tol`.
+* `tol::Float64=1e-15` : ALS converges when fidelity change between two iterations is smaller than `tol`.
 * `check_interval::Int=0` : Set number of iterations to print information. Output is suppressed when `check_interval <= 0`. 
 """
 @kwdef struct ALSTruncation

From 0545de0e81e444284ea8190857ca68e7ff9b8275 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Thu, 13 Nov 2025 15:42:59 +0800
Subject: [PATCH 30/36] Change `gate_bothsides` to `purified`

---
 src/algorithms/time_evolution/evoltools.jl         |  6 ++++++
 src/algorithms/time_evolution/simpleupdate.jl      | 10 +++++++---
 src/algorithms/time_evolution/simpleupdate3site.jl |  6 +++---
 src/algorithms/time_evolution/time_evolve.jl       |  4 ++--
 test/examples/j1j2_finiteT.jl                      |  6 +++---
 test/examples/tf_ising_finiteT.jl                  |  4 ++--
 6 files changed, 23 insertions(+), 13 deletions(-)

diff --git a/src/algorithms/time_evolution/evoltools.jl b/src/algorithms/time_evolution/evoltools.jl
index 9659bc95..e2bd0088 100644
--- a/src/algorithms/time_evolution/evoltools.jl
+++ b/src/algorithms/time_evolution/evoltools.jl
@@ -1,8 +1,14 @@
 const InfiniteState = Union{InfinitePEPS, InfinitePEPO}
 
+"""
+Process the Trotter time step `dt` according to the intended usage.
+"""
 function _get_dt(
         state::InfiniteState, dt::Number, imaginary_time::Bool
     )
+    # PEPS update: exp(-H dt)|ψ⟩
+    # PEPO update (purified): exp(-H dt/2)|ρ⟩
+    # PEPO update (not purified): exp(-H dt/2) ρ exp(-H dt/2)
     dt′ = (state isa InfinitePEPS) ? dt : (dt / 2)
     if (state isa InfinitePEPO)
         @assert size(state)[3] == 1
diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index abb54ea1..26869bf0 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -16,8 +16,12 @@ $(TYPEDFIELDS)
     force_3site::Bool = false
     "(Only applicable to InfinitePEPS) When true, assume bipartite unit cell structure"
     bipartite::Bool = false
-    "(Only applicable to InfinitePEPO) When true (or false), the PEPO is updated as `exp(-H dt/2) * ρ * exp(-H dt/2)` (or `exp(-H dt) ρ`) in each Trotter step"
-    gate_bothsides::Bool = false
+    "(Only applicable to InfinitePEPO) 
+    When true, the PEPO is regarded as a purified PEPS, and updated as
+    `|ρ(t + dt)⟩ = exp(-H dt/2) |ρ(t)⟩`.
+    When false, the PEPO is updated as 
+    `ρ(t + dt) = exp(-H dt/2) ρ(t) exp(-H dt/2)`."
+    purified::Bool = true
 end
 
 """
@@ -158,7 +162,7 @@ function su_iter(
     )
     Nr, Nc, = size(state)
     state2, env2, ϵ = deepcopy(state), deepcopy(env), 0.0
-    gate_axs = alg.gate_bothsides ? (1:2) : (1:1)
+    gate_axs = alg.purified ? (1:1) : (1:2)
     for r in 1:Nr, c in 1:Nc
         term = get_gateterm(gate, (CartesianIndex(r, c), CartesianIndex(r, c + 1)))
         trunc = truncation_strategy(alg.trunc, 1, r, c)
diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl
index 9be45a7f..7b80acde 100644
--- a/src/algorithms/time_evolution/simpleupdate3site.jl
+++ b/src/algorithms/time_evolution/simpleupdate3site.jl
@@ -462,7 +462,7 @@ end
 function _su3site_se!(
         state::InfiniteState, gs::Vector{T}, env::SUWeight,
         row::Int, col::Int, truncs::Vector{E};
-        gate_bothsides::Bool = true
+        purified::Bool = true
     ) where {T <: AbstractTensorMap, E <: TruncationStrategy}
     Nr, Nc = size(state)
     @assert 1 <= row <= Nr && 1 <= col <= Nc
@@ -477,7 +477,7 @@ function _su3site_se!(
     # weights in the cluster
     wt_idxs = ((1, row, col), (2, row, cp1))
     # apply gate MPOs
-    gate_axs = gate_bothsides ? (1:2) : (1:1)
+    gate_axs = purified ? (1:1) : (1:2)
     ϵs = nothing
     for gate_ax in gate_axs
         _apply_gatempo!(Ms, gs; gate_ax)
@@ -525,7 +525,7 @@ function su_iter(
                 truncation_strategy(trunc, 1, r, c)
                 truncation_strategy(trunc, 2, r, _next(c, Nc))
             ]
-            ϵ = _su3site_se!(state2, gs, env2, r, c, truncs; alg.gate_bothsides)
+            ϵ = _su3site_se!(state2, gs, env2, r, c, truncs; alg.purified)
         end
         state2, env2 = rotl90(state2), rotl90(env2)
         trunc = rotl90(trunc)
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index 5071bc1a..9ec72947 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -35,8 +35,8 @@ function _timeevol_sanity_check(
     Nr, Nc, = size(ψ₀)
     @assert (Nr >= 2 && Nc >= 2) "Unit cell size for simple update should be no smaller than (2, 2)."
     @assert Pspaces == physicalspace(ψ₀) "Physical spaces of `ψ₀` do not match `Pspaces`."
-    if hasfield(typeof(alg), :gate_bothsides) && alg.gate_bothsides
-        @assert ψ₀ isa InfinitePEPO "alg.gate_bothsides = true is only compatible with PEPO."
+    if hasfield(typeof(alg), :purified) && !alg.purified
+        @assert ψ₀ isa InfinitePEPO "alg.purified = false is only applicable to PEPO."
     end
     if hasfield(typeof(alg), :bipartite) && alg.bipartite
         @assert Nr == Nc == 2 "`bipartite = true` requires 2 x 2 unit cell size."
diff --git a/test/examples/j1j2_finiteT.jl b/test/examples/j1j2_finiteT.jl
index 1ec19bfb..707d3605 100644
--- a/test/examples/j1j2_finiteT.jl
+++ b/test/examples/j1j2_finiteT.jl
@@ -32,7 +32,7 @@ check_interval = 100
 
 # PEPO approach
 dt, nstep = 1.0e-3, 600
-alg = SimpleUpdate(; trunc = trunc_pepo, gate_bothsides = true)
+alg = SimpleUpdate(; trunc = trunc_pepo, purified = false)
 evolver = TimeEvolver(pepo0, ham, dt, nstep, alg, wts0)
 pepo, wts, info = time_evolve(evolver; check_interval)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
@@ -42,9 +42,9 @@ energy = expectation_value(pepo, ham, env) / (Nr * Nc)
 @test energy ≈ bm[2] atol = 5.0e-3
 
 # PEPS (purified PEPO) approach
-alg = SimpleUpdate(; trunc = trunc_pepo, gate_bothsides = false)
+alg = SimpleUpdate(; trunc = trunc_pepo, purified = true)
 evolver = TimeEvolver(pepo0, ham, dt, nstep, alg, wts0)
-pepo, wts, = time_evolve(evolver; check_interval)
+pepo, wts, info = time_evolve(evolver; check_interval)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 energy = expectation_value(pepo, ham, env) / (Nr * Nc)
 @info "β = $(dt * nstep) / 2: tr(ρH) = $(energy)"
diff --git a/test/examples/tf_ising_finiteT.jl b/test/examples/tf_ising_finiteT.jl
index 109b4109..bb2185c5 100644
--- a/test/examples/tf_ising_finiteT.jl
+++ b/test/examples/tf_ising_finiteT.jl
@@ -60,7 +60,7 @@ dt, nstep = 1.0e-3, 400
 # when g = 2, β = 0.4 and 2β = 0.8 belong to two phases (without and with nonzero σᶻ)
 
 # PEPO approach: results at β, or T = 2.5
-alg = SimpleUpdate(; trunc = trunc_pepo, gate_bothsides = true)
+alg = SimpleUpdate(; trunc = trunc_pepo, purified = false)
 pepo, wts, info = time_evolve(pepo0, ham, dt, nstep, alg, wts0)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 result_β = measure_mag(pepo, env)
@@ -77,7 +77,7 @@ result_2β = measure_mag(pepo, env)
 @test isapprox(abs.(result_2β), bm_2β, rtol = 1.0e-4)
 
 # Purification approach: results at 2β, or T = 1.25
-alg = SimpleUpdate(; trunc = trunc_pepo, gate_bothsides = false)
+alg = SimpleUpdate(; trunc = trunc_pepo, purified = true)
 pepo, wts, info = time_evolve(pepo0, ham, dt, 2 * nstep, alg, wts0)
 env = converge_env(InfinitePEPS(pepo), 8)
 result_2β′ = measure_mag(pepo, env; purified = true)

From c3de149261aa7eb34bb8c15304a9d2856e28c72c Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Thu, 13 Nov 2025 15:44:35 +0800
Subject: [PATCH 31/36] Remove `rand` from finite-T SU tests

---
 test/examples/j1j2_finiteT.jl     | 5 +----
 test/examples/tf_ising_finiteT.jl | 3 ---
 2 files changed, 1 insertion(+), 7 deletions(-)

diff --git a/test/examples/j1j2_finiteT.jl b/test/examples/j1j2_finiteT.jl
index 707d3605..991521c3 100644
--- a/test/examples/j1j2_finiteT.jl
+++ b/test/examples/j1j2_finiteT.jl
@@ -1,12 +1,9 @@
 using Test
-using Random
 using LinearAlgebra
 using TensorKit
 import MPSKitModels: σˣ, σᶻ
 using PEPSKit
 
-Random.seed!(10235876)
-
 # Benchmark energy from high-temperature expansion
 # at β = 0.3, 0.6
 # Physical Review B 86, 045139 (2012) Fig. 15-16
@@ -14,7 +11,7 @@ bm = [-0.1235, -0.213]
 
 function converge_env(state, χ::Int)
     trunc1 = truncrank(χ) & truncerror(; atol = 1.0e-12)
-    env0 = CTMRGEnv(randn, Float64, state, Vect[SU2Irrep](0 => 1))
+    env0 = CTMRGEnv(ones, Float64, state, Vect[SU2Irrep](0 => 1))
     env, = leading_boundary(env0, state; alg = :sequential, trunc = trunc1, tol = 1.0e-10)
     return env
 end
diff --git a/test/examples/tf_ising_finiteT.jl b/test/examples/tf_ising_finiteT.jl
index bb2185c5..37176cc5 100644
--- a/test/examples/tf_ising_finiteT.jl
+++ b/test/examples/tf_ising_finiteT.jl
@@ -1,12 +1,9 @@
 using Test
-using Random
 using LinearAlgebra
 using TensorKit
 import MPSKitModels: σˣ, σᶻ
 using PEPSKit
 
-Random.seed!(10235876)
-
 # Benchmark data of [σx, σz] from HOTRG
 # Physical Review B 86, 045139 (2012) Fig. 15-16
 bm_β = [0.5632, 0.0]

From 7e7a60d747efe10dfcce26680d7feccdee9c49aa Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Fri, 14 Nov 2025 16:05:48 +0800
Subject: [PATCH 32/36] Move `InfiniteState`

---
 src/algorithms/time_evolution/evoltools.jl | 2 --
 src/operators/infinitepepo.jl              | 2 ++
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/algorithms/time_evolution/evoltools.jl b/src/algorithms/time_evolution/evoltools.jl
index e2bd0088..627e38bf 100644
--- a/src/algorithms/time_evolution/evoltools.jl
+++ b/src/algorithms/time_evolution/evoltools.jl
@@ -1,5 +1,3 @@
-const InfiniteState = Union{InfinitePEPS, InfinitePEPO}
-
 """
 Process the Trotter time step `dt` according to the intended usage.
 """
diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl
index e16b4cfc..c1aaa1c3 100644
--- a/src/operators/infinitepepo.jl
+++ b/src/operators/infinitepepo.jl
@@ -27,6 +27,8 @@ struct InfinitePEPO{T <: PEPOTensor}
     end
 end
 
+const InfiniteState = Union{InfinitePEPS, InfinitePEPO}
+
 ## Constructors
 
 """

From aeb9ada29680ccb34bbba1f4c03d819e23f79313 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sat, 15 Nov 2025 10:18:32 +0800
Subject: [PATCH 33/36] Remove functionality to reset `iter` from `timestep`

---
 src/algorithms/time_evolution/simpleupdate.jl | 11 +++--------
 test/timeevol/timestep.jl                     |  2 +-
 2 files changed, 4 insertions(+), 9 deletions(-)

diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index 26869bf0..f34428e8 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -206,22 +206,17 @@ end
 
 """
     timestep(
-        it::TimeEvolver{<:SimpleUpdate}, ψ::InfiniteState, env::SUWeight;
-        iter::Int = it.state.iter, t::Float64 = it.state.t
+        it::TimeEvolver{<:SimpleUpdate}, ψ::InfiniteState, env::SUWeight
     ) -> (ψ, env, info)
 
 Given the TimeEvolver iterator `it`, perform one step of time evolution
 on the input state `ψ` and its environment `env`.
-
-- Using `iter` and `t` to reset the current iteration number and evolved time
-    respectively of the TimeEvolver `it`.
 """
 function MPSKit.timestep(
-        it::TimeEvolver{<:SimpleUpdate}, ψ::InfiniteState, env::SUWeight;
-        iter::Int = it.state.iter, t::Float64 = it.state.t
+        it::TimeEvolver{<:SimpleUpdate}, ψ::InfiniteState, env::SUWeight
     )
     _timeevol_sanity_check(ψ, physicalspace(it.state.ψ), it.alg)
-    state = SUState(iter, t, ψ, env)
+    state = SUState(it.state.iter, it.state.t, ψ, env)
     result = iterate(it, state)
     if result === nothing
         @warn "TimeEvolver `it` has already reached the end."
diff --git a/test/timeevol/timestep.jl b/test/timeevol/timestep.jl
index c1513bd9..39153a89 100644
--- a/test/timeevol/timestep.jl
+++ b/test/timeevol/timestep.jl
@@ -15,7 +15,7 @@ using PEPSKit
     evolver = TimeEvolver(ψ0, H, dt, nstep, alg, env0)
     ψ1, env1, info1 = deepcopy(ψ0), deepcopy(env0), nothing
     for iter in 0:(nstep - 1)
-        ψ1, env1, info1 = timestep(evolver, ψ1, env1; iter)
+        ψ1, env1, info1 = timestep(evolver, ψ1, env1)
     end
     # time_evolve
     ψ2, env2, info2 = time_evolve(ψ0, H, dt, nstep, alg, env0)

From 2e3c6c2ad1f881abcd15147af498ca80ec8b2ff8 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sat, 15 Nov 2025 10:27:01 +0800
Subject: [PATCH 34/36] Docstring updates

---
 examples/hubbard_su/main.jl                  | 4 ++--
 src/algorithms/time_evolution/time_evolve.jl | 3 ++-
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/examples/hubbard_su/main.jl b/examples/hubbard_su/main.jl
index 77f13338..786dab68 100644
--- a/examples/hubbard_su/main.jl
+++ b/examples/hubbard_su/main.jl
@@ -76,8 +76,8 @@ md"""
 In order to compute the energy expectation value with evolved PEPS, we need to converge a
 CTMRG environment on it. We first converge an environment with a small enviroment dimension,
 which is initialized using the simple update bond weights. Next we use it to initialize 
-another run with bigger environment dimension. The dynamic adjustion of environment dimension
-is achieved by using `trunc=truncrank(χ)` in the CTMRG runs:
+another run with bigger environment dimension. The dynamic adjustment of environment dimension
+is achieved by using `trunc=truncrank(χ)` with different `χ`s in the CTMRG runs:
 """
 
 χenv₀, χenv = 6, 16
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index 9ec72947..807caaa0 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -23,7 +23,8 @@ mutable struct TimeEvolver{TE <: TimeEvolution, G, S, N <: Number}
     nstep::Int
     "Trotter gates"
     gate::G
-    "PEPS/PEPO and its environment"
+    "Internal state of the iterator, including the number of
+    already performed iterations, evolved time, PEPS/PEPO and its environment"
     state::S
 end
 

From 45538386b0329e2fef788d6309d8bcb45b291dfb Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sat, 15 Nov 2025 10:36:39 +0800
Subject: [PATCH 35/36] Reduce unicode usage

---
 src/algorithms/time_evolution/simpleupdate.jl | 68 +++++++++----------
 1 file changed, 34 insertions(+), 34 deletions(-)

diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index f34428e8..17e3fc46 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -33,35 +33,35 @@ struct SUState{S <: InfiniteState, E <: SUWeight, N <: Number}
     "evolved time"
     t::N
     "PEPS/PEPO"
-    ψ::S
+    psi::S
     "SUWeight environment"
     env::E
 end
 
 """
     TimeEvolver(
-        ψ₀::InfiniteState, H::LocalOperator, dt::Number, nstep::Int, 
-        alg::SimpleUpdate, env₀::SUWeight; t₀::Number = 0.0
+        psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int, 
+        alg::SimpleUpdate, env0::SUWeight; t0::Number = 0.0
     )
 
 Initialize a TimeEvolver with Hamiltonian `H` and simple update `alg`, 
-starting from the initial state `ψ₀` and SUWeight environment `env₀`.
+starting from the initial state `psi0` and SUWeight environment `env0`.
 
-- The initial time is specified by `t₀`.
+- The initial time is specified by `t0`.
 """
 function TimeEvolver(
-        ψ₀::InfiniteState, H::LocalOperator, dt::Number, nstep::Int,
-        alg::SimpleUpdate, env₀::SUWeight; t₀::Number = 0.0
+        psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int,
+        alg::SimpleUpdate, env0::SUWeight; t0::Number = 0.0
     )
     # sanity checks
-    _timeevol_sanity_check(ψ₀, physicalspace(H), alg)
+    _timeevol_sanity_check(psi0, physicalspace(H), alg)
     # create Trotter gates
     nnonly = is_nearest_neighbour(H)
     use_3site = alg.force_3site || !nnonly
     if alg.bipartite
         @assert !use_3site "3-site simple update is incompatible with bipartite lattice."
     end
-    dt′ = _get_dt(ψ₀, dt, alg.imaginary_time)
+    dt′ = _get_dt(psi0, dt, alg.imaginary_time)
     gate = if use_3site
         [
             _get_gatempos_se(H, dt′),
@@ -72,7 +72,7 @@ function TimeEvolver(
     else
         get_expham(H, dt′)
     end
-    state = SUState(0, t₀, ψ₀, env₀)
+    state = SUState(0, t0, psi0, env0)
     return TimeEvolver(alg, dt, nstep, gate, state)
 end
 
@@ -195,28 +195,28 @@ end
 function Base.iterate(it::TimeEvolver{<:SimpleUpdate}, state = it.state)
     iter, t = state.iter, state.t
     (iter == it.nstep) && return nothing
-    ψ, env, ϵ = su_iter(state.ψ, it.gate, it.alg, state.env)
+    psi, env, ϵ = su_iter(state.psi, it.gate, it.alg, state.env)
     # update internal state
     iter += 1
     t += it.dt
-    it.state = SUState(iter, t, ψ, env)
+    it.state = SUState(iter, t, psi, env)
     info = (; t, ϵ)
-    return (ψ, env, info), it.state
+    return (psi, env, info), it.state
 end
 
 """
     timestep(
-        it::TimeEvolver{<:SimpleUpdate}, ψ::InfiniteState, env::SUWeight
-    ) -> (ψ, env, info)
+        it::TimeEvolver{<:SimpleUpdate}, psi::InfiniteState, env::SUWeight
+    ) -> (psi, env, info)
 
 Given the TimeEvolver iterator `it`, perform one step of time evolution
-on the input state `ψ` and its environment `env`.
+on the input state `psi` and its environment `env`.
 """
 function MPSKit.timestep(
-        it::TimeEvolver{<:SimpleUpdate}, ψ::InfiniteState, env::SUWeight
+        it::TimeEvolver{<:SimpleUpdate}, psi::InfiniteState, env::SUWeight
     )
-    _timeevol_sanity_check(ψ, physicalspace(it.state.ψ), it.alg)
-    state = SUState(it.state.iter, it.state.t, ψ, env)
+    _timeevol_sanity_check(psi, physicalspace(it.state.psi), it.alg)
+    state = SUState(it.state.iter, it.state.t, psi, env)
     result = iterate(it, state)
     if result === nothing
         @warn "TimeEvolver `it` has already reached the end."
@@ -230,7 +230,7 @@ end
     time_evolve(
         it::TimeEvolver{<:SimpleUpdate}; 
         tol::Float64 = 0.0, check_interval::Int = 500
-    ) -> (ψ, env, info)
+    ) -> (psi, env, info)
 
 Perform time evolution to the end of TimeEvolver iterator `it`,
 or until convergence of SUWeight set by a positive `tol`.
@@ -246,10 +246,10 @@ function MPSKit.time_evolve(
     time_start = time()
     check_convergence = (tol > 0)
     if check_convergence
-        @assert (it.state.ψ isa InfinitePEPS) && it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking."
+        @assert (it.state.psi isa InfinitePEPS) && it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking."
     end
     env0, time0 = it.state.env, time()
-    for (ψ, env, info) in it
+    for (psi, env, info) in it
         iter = it.state.iter
         diff = compare_weights(env0, env)
         stop = (iter == it.nstep) || (diff < tol)
@@ -275,7 +275,7 @@ function MPSKit.time_evolve(
         if stop
             time_end = time()
             @info @sprintf("Simple update finished. Total time elasped: %.2f s", time_end - time_start)
-            return ψ, env, info
+            return psi, env, info
         else
             env0 = env
         end
@@ -286,27 +286,27 @@ end
 
 """
     time_evolve(
-        ψ₀::Union{InfinitePEPS, InfinitePEPO}, H::LocalOperator, 
-        dt::Number, nstep::Int, alg::SimpleUpdate, env₀::SUWeight;
-        tol::Float64 = 0.0, t₀::Number = 0.0, check_interval::Int = 500
-    ) -> (ψ, env, info)
+        psi0::Union{InfinitePEPS, InfinitePEPO}, H::LocalOperator, 
+        dt::Number, nstep::Int, alg::SimpleUpdate, env0::SUWeight;
+        tol::Float64 = 0.0, t0::Number = 0.0, check_interval::Int = 500
+    ) -> (psi, env, info)
 
-Perform time evolution on the initial state `ψ₀` and initial environment `env₀`
+Perform time evolution on the initial state `psi0` and initial environment `env0`
 with Hamiltonian `H`, using SimpleUpdate algorithm `alg`, time step `dt` for 
 `nstep` number of steps. 
 
 - Setting `tol > 0` enables convergence check (for imaginary time evolution of InfinitePEPS only).
     For other usages it should not be changed.
-- Use `t₀` to specify the initial time of `ψ₀`.
+- Use `t0` to specify the initial time of the evolution.
 - `check_interval` sets the interval to output information. Output during the evolution can be turned off by setting `check_interval <= 0`.
 - `info` is a NamedTuple containing information of the evolution, 
-    including the time `info.t` evolved since `ψ₀`.
+    including the time `info.t` evolved since `psi0`.
 """
 function MPSKit.time_evolve(
-        ψ₀::InfiniteState, H::LocalOperator, dt::Number, nstep::Int,
-        alg::SimpleUpdate, env₀::SUWeight;
-        tol::Float64 = 0.0, t₀::Number = 0.0, check_interval::Int = 500
+        psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int,
+        alg::SimpleUpdate, env0::SUWeight;
+        tol::Float64 = 0.0, t0::Number = 0.0, check_interval::Int = 500
     )
-    it = TimeEvolver(ψ₀, H, dt, nstep, alg, env₀; t₀)
+    it = TimeEvolver(psi0, H, dt, nstep, alg, env0; t0)
     return time_evolve(it; tol, check_interval)
 end

From 700378f0830472029399553e1bebc3f497d72f9d Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sat, 15 Nov 2025 11:23:58 +0800
Subject: [PATCH 36/36] Fix test

---
 test/examples/tf_ising_finiteT.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/examples/tf_ising_finiteT.jl b/test/examples/tf_ising_finiteT.jl
index 37176cc5..584b21bf 100644
--- a/test/examples/tf_ising_finiteT.jl
+++ b/test/examples/tf_ising_finiteT.jl
@@ -66,7 +66,7 @@ result_β = measure_mag(pepo, env)
 @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2)
 
 # continue to get results at 2β, or T = 1.25
-pepo, wts, info = time_evolve(pepo, ham, dt, nstep, alg, wts; t₀ = β)
+pepo, wts, info = time_evolve(pepo, ham, dt, nstep, alg, wts; t0 = β)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 result_2β = measure_mag(pepo, env)
 @info "tr(σ(x,z)ρ) at T = $(1 / (2β))" result_2β