diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 3c950af33..4d2f583a1 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -105,9 +105,9 @@ include("algorithms/truncation/truncationschemes.jl")
 include("algorithms/truncation/fullenv_truncation.jl")
 include("algorithms/truncation/bond_truncation.jl")
 
-include("algorithms/time_evolution/trotter_gate.jl")
 include("algorithms/time_evolution/apply_gate.jl")
 include("algorithms/time_evolution/apply_mpo.jl")
+include("algorithms/time_evolution/trotter_gate.jl")
 include("algorithms/time_evolution/time_evolve.jl")
 include("algorithms/time_evolution/simpleupdate.jl")
 include("algorithms/time_evolution/simpleupdate3site.jl")
diff --git a/src/algorithms/bp/beliefpropagation.jl b/src/algorithms/bp/beliefpropagation.jl
index 0246d7e20..0aa900912 100644
--- a/src/algorithms/bp/beliefpropagation.jl
+++ b/src/algorithms/bp/beliefpropagation.jl
@@ -59,7 +59,7 @@ function leading_boundary(env₀::BPEnv, network::InfiniteSquareNetwork, alg::Be
 end
 function leading_boundary(env₀::BPEnv, state::InfiniteState, alg::BeliefPropagation)
     if alg.bipartite
-        @assert _state_bipartite_check(state)
+        _is_bipartite(state) || error("Input state is not bipartite")
     end
     return leading_boundary(env₀, InfiniteSquareNetwork(state), alg)
 end
diff --git a/src/algorithms/bp/gaugefix.jl b/src/algorithms/bp/gaugefix.jl
index 561dd2183..56e5e9f30 100644
--- a/src/algorithms/bp/gaugefix.jl
+++ b/src/algorithms/bp/gaugefix.jl
@@ -7,16 +7,6 @@ Algorithm for gauging PEPS with belief propagation fixed point messages.
     # TODO: add options
 end
 
-function _bpenv_bipartite_check(env::BPEnv)
-    for (r, c) in Iterators.product(1:2, 1:2)
-        r′, c′ = _next(r, 2), _next(c, 2)
-        if !all(env[:, r, c] .== env[:, r′, c′])
-            return false
-        end
-    end
-    return true
-end
-
 """
     gauge_fix(psi::Union{InfinitePEPS, InfinitePEPO}, alg::BPGauge, env::BPEnv)
 
@@ -25,7 +15,7 @@ an [`InfinitePEPO`](@ref) interpreted as purified state with two physical legs)
 using fixed point environment `env` of belief propagation.
 """
 function gauge_fix(psi::InfinitePEPS, alg::BPGauge, env::BPEnv)
-    bipartite = _state_bipartite_check(psi) && _bpenv_bipartite_check(env)
+    bipartite = _is_bipartite(psi) && _is_bipartite(env)
     psi′ = copy(psi)
     XXinv = map(eachcoordinate(psi, 1:2)) do I
         _, X, Xinv = _bp_gauge_fix!(CartesianIndex(I), psi′, env)
@@ -52,7 +42,7 @@ function gauge_fix(psi::InfinitePEPO, alg::BPGauge, env::BPEnv)
     Fs = map(Base.Fix2(getindex, 2), psi_Fs)
     psi′, XXinv = gauge_fix(InfinitePEPS(psi′), alg, env)
     # convert back to iPEPO
-    psi′ = map(zip(psi′.A, Fs)) do (t, F)
+    psi′ = map(psi′.A, Fs) do t, F
         return F' * t
     end
     psi′ = reshape(psi′, (Nr, Nc, 1))
diff --git a/src/algorithms/time_evolution/apply_gate.jl b/src/algorithms/time_evolution/apply_gate.jl
index 27f9da61f..0bc331895 100644
--- a/src/algorithms/time_evolution/apply_gate.jl
+++ b/src/algorithms/time_evolution/apply_gate.jl
@@ -1,3 +1,5 @@
+const NNGate{T, S} = AbstractTensorMap{T, S, 2, 2}
+
 """
 Apply 1-site `gate` on the PEPS or PEPO tensor `a`.
 """
@@ -125,8 +127,8 @@ Apply 2-site `gate` on the reduced matrices `a`, `b`
 """
 function _apply_gate(
         a::AbstractTensorMap, b::AbstractTensorMap,
-        gate::AbstractTensorMap{T, S, 2, 2}, trunc::TruncationStrategy
-    ) where {T <: Number, S <: ElementarySpace}
+        gate::NNGate, trunc::TruncationStrategy
+    )
     V = space(b, 1)
     need_flip = isdual(V)
     if isdual(space(a, 2))
@@ -134,11 +136,6 @@ function _apply_gate(
     else
         @tensor a2b2[-1 -2; -3 -4] := gate[-2 -3; 1 2] * a[-1 1 3] * b[3 2 -4]
     end
-    trunc = if trunc isa FixedSpaceTruncation
-        need_flip ? truncspace(flip(V)) : truncspace(V)
-    else
-        trunc
-    end
     a, s, b, ϵ = svd_trunc!(a2b2; trunc)
     a, b = absorb_s(a, s, b)
     if need_flip
diff --git a/src/algorithms/time_evolution/apply_mpo.jl b/src/algorithms/time_evolution/apply_mpo.jl
index b7ae6f42f..bb11b805c 100644
--- a/src/algorithms/time_evolution/apply_mpo.jl
+++ b/src/algorithms/time_evolution/apply_mpo.jl
@@ -229,14 +229,8 @@ function _get_allprojs(
     N = length(Ms)
     Rs, Ls = _get_allRLs(Ms)
     @assert length(truncs) == N - 1
-    projs_errs = map(1:(N - 1)) do i
-        trunc = if isa(truncs[i], FixedSpaceTruncation)
-            tspace = space(Ms[i + 1], 1)
-            isdual(tspace) ? truncspace(flip(tspace)) : truncspace(tspace)
-        else
-            truncs[i]
-        end
-        return _proj_from_RL(Rs[i], Ls[i]; trunc)
+    projs_errs = map(Rs, Ls, truncs) do R, L, trunc
+        return _proj_from_RL(R, L; trunc)
     end
     Pas = map(Base.Fix2(getindex, 1), projs_errs)
     wts = map(Base.Fix2(getindex, 2), projs_errs)
diff --git a/src/algorithms/time_evolution/gaugefix_su.jl b/src/algorithms/time_evolution/gaugefix_su.jl
index 324a8d604..13a2782d6 100644
--- a/src/algorithms/time_evolution/gaugefix_su.jl
+++ b/src/algorithms/time_evolution/gaugefix_su.jl
@@ -37,22 +37,24 @@ end
 Fix the gauge of `psi` using trivial simple update.
 """
 function gauge_fix(psi::InfiniteState, alg::SUGauge)
+    time0 = time()
     gates = _trivial_gates(scalartype(psi), physicalspace(psi))
-    su_alg = SimpleUpdate(; trunc = FixedSpaceTruncation(), bipartite = _state_bipartite_check(psi))
+    trunc = _get_fixedspacetrunc(psi)
+    su_alg = SimpleUpdate(; trunc, bipartite = _is_bipartite(psi))
     wts0 = SUWeight(psi)
     # use default constructor to avoid calculation of exp(-H * 0)
     evolver = TimeEvolver(su_alg, 0.0, alg.maxiter, gates, SUState(0, 0.0, psi, wts0))
     for (i, (psi′, wts, info)) in enumerate(evolver)
         ϵ = compare_weights(wts, wts0)
         if i >= alg.miniter && ϵ < alg.tol
-            @info "Trivial SU conv $i: |Δλ| = $ϵ."
+            @info "Trivial SU conv $i: |Δλ| = $ϵ, time = $(time() - time0) s"
             return psi′, wts, ϵ
         end
         if i == alg.maxiter
-            @warn "Trivial SU cancel $i: |Δλ| = $ϵ."
+            @warn "Trivial SU cancel $i: |Δλ| = $ϵ, time = $(time() - time0) s"
             return psi′, wts, ϵ
         end
-        wts0 = deepcopy(wts)
+        wts0 = wts
     end
     return
 end
diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index a28f86d8b..75bd72a2e 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -7,19 +7,19 @@ Algorithm struct for simple update (SU) of InfinitePEPS or InfinitePEPO.
 
 $(TYPEDFIELDS)
 """
-@kwdef struct SimpleUpdate <: TimeEvolution
+@kwdef struct SimpleUpdate{T <: TruncationStrategy} <: TimeEvolution
     "Truncation strategy for bonds updated by Trotter gates"
-    trunc::TruncationStrategy
+    trunc::T
     "When true (or false), the Trotter gate is `exp(-H dt)` (or `exp(-iH dt)`)"
     imaginary_time::Bool = true
     "When true, force decomposition of nearest neighbor gates to MPOs."
     force_mpo::Bool = false
     "When true, assume bipartite unit cell structure"
     bipartite::Bool = false
-    "(Only applicable to InfinitePEPO) 
+    "(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 
+    When false, the PEPO is updated as
     `ρ(t + dt) = exp(-H dt/2) ρ(t) exp(-H dt/2)`."
     purified::Bool = true
 end
@@ -62,6 +62,12 @@ function TimeEvolver(
     # create Trotter gates
     gate = trotterize(H, dt′; symmetrize_gates, force_mpo = alg.force_mpo)
     state = SUState(0, t0, psi0, env0)
+    # convert FixedSpaceTruncation to site-dependent `truncspace`s
+    if alg.trunc isa FixedSpaceTruncation
+        trunc = _get_fixedspacetrunc(psi0)
+        @reset alg.trunc = trunc
+    end
+    # TODO: bipartite check for alg.trunc after equality is defined for all kinds of truncation strategies
     # TODO: check gates for bipartite case
     return TimeEvolver(alg, dt, nstep, gate, state)
 end
@@ -89,8 +95,8 @@ function _su_iter!(
         sites::Vector{CartesianIndex{2}}, alg::SimpleUpdate
     )
     Nr, Nc = size(state)
-    truncs = _get_cluster_trunc(alg.trunc, sites, (Nr, Nc))
-    @assert length(sites) == 2 && length(truncs) == 1
+    @assert length(sites) == 2
+    trunc = only(_get_cluster_trunc(alg.trunc, sites, (Nr, Nc)))
     Ms, open_vaxs, = _get_cluster(state, sites, env; permute = false)
     normalize!.(Ms, Inf)
     # rotate
@@ -101,7 +107,7 @@ function _su_iter!(
     gate_axs = alg.purified ? (1:1) : (1:2)
     for gate_ax in gate_axs
         X, a, b, Y = _qr_bond(A, B; gate_ax, positive = true)
-        a, s, b, ϵ′ = _apply_gate(a, b, gate, truncs[1])
+        a, s, b, ϵ′ = _apply_gate(a, b, gate, trunc)
         ϵ = max(ϵ, ϵ′)
         A, B = _qr_bond_undo(X, a, b, Y)
     end
@@ -148,14 +154,14 @@ function su_iter(
             (!alg.bipartite) && continue
             if d == 1
                 rp1, cp1 = _next(r, Nr), _next(c, Nc)
-                state2[rp1, cp1] = deepcopy(state2[r, c])
-                state2[rp1, c] = deepcopy(state2[r, cp1])
-                env2[1, rp1, cp1] = deepcopy(env2[1, r, c])
+                state2[rp1, cp1] = copy(state2[r, c])
+                state2[rp1, c] = copy(state2[r, cp1])
+                env2[1, rp1, cp1] = copy(env2[1, r, c])
             else
                 rm1, cm1 = _prev(r, Nr), _prev(c, Nc)
-                state2[rm1, cm1] = deepcopy(state2[r, c])
-                state2[r, cm1] = deepcopy(state2[rm1, c])
-                env2[2, rm1, cm1] = deepcopy(env2[2, r, c])
+                state2[rm1, cm1] = copy(state2[r, c])
+                state2[r, cm1] = copy(state2[rm1, c])
+                env2[2, rm1, cm1] = copy(env2[2, r, c])
             end
         else
             # N-site MPO gate (N ≥ 2)
@@ -202,10 +208,8 @@ function MPSKit.timestep(
 end
 
 """
-    time_evolve(
-        it::TimeEvolver{<:SimpleUpdate}; 
-        tol::Float64 = 0.0, check_interval::Int = 500
-    ) -> (psi, env, info)
+    time_evolve(it; check_interval = 500) -> (psi, env, info)
+    time_evolve(it, H; tol = 1.0e-8, check_interval = 500) -> (psi, env, info)
 
 Perform time evolution to the end of `TimeEvolver` iterator `it`,
 or until convergence of `SUWeight` set by a positive `tol`.
@@ -215,15 +219,41 @@ or until convergence of `SUWeight` set by a positive `tol`.
 - `check_interval` sets the number of iterations between outputs of information.
 """
 function MPSKit.time_evolve(
-        it::TimeEvolver{<:SimpleUpdate};
-        tol::Float64 = 0.0, check_interval::Int = 500
+        it::TimeEvolver{<:SimpleUpdate}; check_interval::Int = 500
     )
     time_start = time()
-    check_convergence = (tol > 0)
     @info "--- Time evolution (simple update), dt = $(it.dt) ---"
-    if check_convergence
-        @assert (it.state.psi isa InfinitePEPS) && it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking."
+    env0, time0 = it.state.env, time()
+    for (psi, env, info) in it
+        iter = it.state.iter
+        diff = compare_weights(env0, env)
+        stop = (iter == it.nstep)
+        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))"
+            @info @sprintf("SU iter %-7d: |Δλ| = %.3e. Time = %.3f s/it", iter, diff, time1 - time0)
+        end
+        if stop
+            time_end = time()
+            @info @sprintf("Time evolution finished in %.2f s", time_end - time_start)
+            return psi, env, info
+        else
+            env0 = env
+        end
+        time0 = time()
     end
+    return
+end
+
+function MPSKit.time_evolve(
+        it::TimeEvolver{<:SimpleUpdate, G, S}, H::LocalOperator;
+        tol::Float64 = 1.0e-8, check_interval::Int = 500
+    ) where {G, S <: SUState{<:InfinitePEPS}}
+    time_start = time()
+    @info "--- Time evolution (simple update), dt = $(it.dt) ---"
+    @assert it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking."
     env0, time0 = it.state.env, time()
     for (psi, env, info) in it
         iter = it.state.iter
@@ -233,16 +263,20 @@ function MPSKit.time_evolve(
             ((iter % check_interval == 0) || (iter == 1) || stop)
         time1 = time()
         if showinfo
+            # TODO: convert to BPEnv instead
+            ctmenv = CTMRGEnv(env)
+            energy = real(expectation_value(psi, H, ctmenv)) / prod(size(psi))
             @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))"
-            @info @sprintf("SU iter %-7d: |Δλ| = %.3e. Time = %.3f s/it", iter, diff, time1 - time0)
+            @info @sprintf(
+                "SU iter %-7d: E ≈ %.5f, |Δλ| = %.3e. Time = %.3f s/it",
+                iter, energy, 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
+        if (iter == it.nstep) && (diff >= tol)
+            @warn "SU: bond weights have not converged."
+        end
+        if diff < tol
+            @info "SU: bond weights have converged."
         end
         if stop
             time_end = time()
@@ -269,7 +303,6 @@ algorithm `alg`, time step `dt` for `nstep` number of steps.
 
 - Set `symmetrize_gates = true` for second-order Trotter decomposition.
 - Set `tol > 0` to enable convergence check (for imaginary time evolution of iPEPS only).
-    For other usages it should not be changed.
 - 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, 
@@ -281,5 +314,9 @@ function MPSKit.time_evolve(
         tol::Float64 = 0.0, t0::Number = 0.0, check_interval::Int = 500
     )
     it = TimeEvolver(psi0, H, dt, nstep, alg, env0; t0, symmetrize_gates)
-    return time_evolve(it; tol, check_interval)
+    return if tol == 0
+        time_evolve(it; check_interval)
+    else
+        time_evolve(it, H; tol, check_interval)
+    end
 end
diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl
index 4b091c830..d306c68db 100644
--- a/src/algorithms/time_evolution/simpleupdate3site.jl
+++ b/src/algorithms/time_evolution/simpleupdate3site.jl
@@ -147,7 +147,7 @@ function _get_cluster(
         p = invperm((p1..., p2...))
         return (p[1:Np], p[(Np + 1):end])
     end
-    Ms = map(zip(sites, open_vaxs, perms)) do (site, vaxs, perm)
+    Ms = map(sites, open_vaxs, perms) do site, vaxs, perm
         s = CartesianIndex(mod1(site[1], Nr), mod1(site[2], Nc))
         M = if env === nothing
             state[s]
@@ -190,7 +190,7 @@ function _su_iter!(
     # restore virtual arrows in `Ms`
     _flip_virtuals!(Ms, flips)
     # update env weights
-    bond_revs = map(zip(sites, Iterators.drop(sites, 1))) do (site1, site2)
+    bond_revs = map(sites, Iterators.drop(sites, 1)) do site1, site2
         _nn_bondrev(site1, site2, (Nr, Nc))
     end
     for (wt, (bond, rev), flip) in zip(wts, bond_revs, flips)
@@ -217,24 +217,14 @@ updated by the Trotter evolution MPO.
 """
 function _get_cluster_trunc(
         trunc::TruncationStrategy, sites::Vector{CartesianIndex{2}},
-        (Nrow, Ncol)::NTuple{2, Int}
+        unitcell::NTuple{2, Int}
     )
-    return map(zip(sites, Iterators.drop(sites, 1))) do (site1, site2)
-        diff = site2 - site1
-        if diff == CartesianIndex(0, 1)
-            r, c = mod1(site1[1], Nrow), mod1(site1[2], Ncol)
-            return truncation_strategy(trunc, 1, r, c)
-        elseif diff == CartesianIndex(0, -1)
-            r, c = mod1(site2[1], Nrow), mod1(site2[2], Ncol)
-            return truncation_strategy(trunc, 1, r, c)
-        elseif diff == CartesianIndex(1, 0)
-            r, c = mod1(site2[1], Nrow), mod1(site2[2], Ncol)
-            return truncation_strategy(trunc, 2, r, c)
-        elseif diff == CartesianIndex(-1, 0)
-            r, c = mod1(site1[1], Nrow), mod1(site1[2], Ncol)
-            return truncation_strategy(trunc, 2, r, c)
-        else
-            error("The path `sites` contains a long-range bond.")
+    return map(sites, Iterators.drop(sites, 1)) do site1, site2
+        (d, r, c), rev = _nn_bondrev(site1, site2, unitcell)
+        t = truncation_strategy(trunc, d, r, c)
+        if rev && isa(t, TruncationSpace)
+            t = truncspace(flip(t.space)')
         end
+        return t
     end
 end
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index ee5f236c1..c1cc50537 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -30,22 +30,6 @@ end
 
 Base.iterate(it::TimeEvolver) = iterate(it, it.state)
 
-function _state_bipartite_check(psi::InfiniteState)
-    if isa(psi, InfinitePEPO)
-        @assert size(psi, 3) == 1 "Input InfinitePEPO is expected to have only one layer."
-    end
-    if !(size(psi, 1) == size(psi, 2) == 2)
-        return false
-    end
-    for (r, c) in Iterators.product(1:2, 1:2)
-        r′, c′ = _next(r, 2), _next(c, 2)
-        if psi[r, c] != psi[r′, c′]
-            return false
-        end
-    end
-    return true
-end
-
 """
 Process the Trotter time step `dt` according to the intended usage.
 """
@@ -79,7 +63,7 @@ function _timeevol_sanity_check(
         @assert ψ₀ isa InfinitePEPO "alg.purified = false is only applicable to PEPO."
     end
     if hasfield(typeof(alg), :bipartite) && alg.bipartite
-        @assert _state_bipartite_check(ψ₀) "Input state is not bipartite with 2 x 2 unit cell."
+        @assert _is_bipartite(ψ₀) "Input state is not bipartite with 2 x 2 unit cell."
     end
     return nothing
 end
@@ -94,3 +78,21 @@ function MPSKit.infinite_temperature_density_matrix(H::LocalOperator)
     end
     return InfinitePEPO(cat(A; dims = 3))
 end
+
+"""
+Get the `SiteDependentTruncation` used by time evolution
+that preserves virtual spaces of `state`.
+"""
+function _get_fixedspacetrunc(state::InfiniteState)
+    if state isa InfinitePEPO
+        size(state, 3) != 1 && error("Input InfinitePEPO is expected to have only one layer.")
+    end
+    Nr, Nc = size(state)
+    return SiteDependentTruncation(
+        map(Iterators.product(1:2, 1:Nr, 1:Nc)) do (d, r, c)
+            V = domain(state[r, c], (d == 1) ? EAST : NORTH)
+            isdual(V) && (V = flip(V))
+            return truncspace(V)
+        end
+    )
+end
diff --git a/src/algorithms/time_evolution/trotter_gate.jl b/src/algorithms/time_evolution/trotter_gate.jl
index 24248ad9c..67ae916d3 100644
--- a/src/algorithms/time_evolution/trotter_gate.jl
+++ b/src/algorithms/time_evolution/trotter_gate.jl
@@ -1,5 +1,3 @@
-const NNGate{T, S} = AbstractTensorMap{T, S, 2, 2}
-
 """
 Convert an N-site gate (N ≥ 2) to MPO by SVD, 
 in which the axes are ordered as
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index 647ad1221..d9b7802a3 100644
--- a/src/algorithms/toolbox.jl
+++ b/src/algorithms/toolbox.jl
@@ -357,7 +357,7 @@ function product_peps(peps_args...; unitcell = (1, 1), noise_amp = 1.0e-2, state
         all(dim.(space.(noise_peps.A, 1)) .== length.(state_vector)) ||
             throw(ArgumentError("state vectors must match the physical dimension"))
     end
-    prod_tensors = map(zip(noise_peps.A, state_vector)) do (t, v)
+    prod_tensors = map(noise_peps.A, state_vector) do t, v
         pt = zero(t)
         pt[][:, 1, 1, 1, 1] .= v
         return pt
diff --git a/src/algorithms/truncation/truncationschemes.jl b/src/algorithms/truncation/truncationschemes.jl
index 71b5e9dbe..7cad1dfd6 100644
--- a/src/algorithms/truncation/truncationschemes.jl
+++ b/src/algorithms/truncation/truncationschemes.jl
@@ -1,16 +1,41 @@
 """
 $(TYPEDEF)
 
-CTMRG specific truncation strategy for `svd_trunc` which keeps the bond space on which the SVD
-is performed fixed. Since different environment directions and unit cell entries might
-have different spaces, this truncation style is different from `TruncationSpace`.
+SVD truncation strategy which preserves the `CTMRGEnv` environment virtual spaces,
+or `InfinitePEPS`, `InfinitePEPO` virtual spaces.
 """
 struct FixedSpaceTruncation <: TruncationStrategy end
 
+"""
+$(TYPEDEF)
+
+SVD truncation strategy specified for each nearest neighbor bond 
+in `InfinitePEPS`, `InfinitePEPO`:
+- `trunc[1, r, c]` applies to the x-bond between `[r, c]` and `[r, c+1]`.
+    If it is a `TruncationSpace`, the space refers to the east domain of
+    `[r, c]` or its `flip`, whichever is non-dual.
+- `trunc[2, r, c]` applies to the y-bond between `[r, c]` and `[r-1, c]`.
+    If it is a `TruncationSpace`, the space refers to the north domain of
+    `[r, c]` or its `flip`, whichever is non-dual.
+"""
 struct SiteDependentTruncation{T <: TruncationStrategy} <: TruncationStrategy
     truncs::Array{T, 3}
+
+    function SiteDependentTruncation(truncs::Array{T, 3}) where {T}
+        # TODO: generalize it to CTMRGEnv
+        size(truncs, 1) != 2 && throw(
+            DimensionMismatch(
+                "The first dimension of `truncs` must have a size of 2. Got $(size(truncs, 1))."
+            )
+        )
+        return new{T}(truncs)
+    end
 end
 
+Base.getindex(trunc::SiteDependentTruncation, args...) = Base.getindex(trunc.truncs, args...)
+
+# TODO: _is_bipartite(trunc::SiteDependentTruncation)
+
 const TRUNCATION_STRATEGY_SYMBOLS = IdDict{Symbol, Type{<:TruncationStrategy}}(
     :notrunc => MatrixAlgebraKit.NoTruncation,
     :truncerror => MatrixAlgebraKit.TruncationByError,
@@ -40,28 +65,68 @@ end
 function truncation_strategy(
         trunc::SiteDependentTruncation, direction::Int, row::Int, col::Int
     )
-    return trunc.truncs[direction, row, col]
+    return trunc[direction, row, col]
 end
 
-# TODO: type piracy
-Base.rotl90(trunc::TruncationStrategy) = trunc
-
-function Base.rotl90(trunc::SiteDependentTruncation)
-    directions, rows, cols = size(trunc.truncs)
-    truncs_rotated = similar(trunc.truncs, directions, cols, rows)
+# rotation of SiteDependentTruncation
+# (similar to rotation of SUWeight)
+function _rotl90_trunc_x(trunc_x::AbstractMatrix{<:TruncationStrategy})
+    trunc_y = rotl90(trunc_x)
+    return trunc_y
+end
+function _rotr90_trunc_x(trunc_x::AbstractMatrix{<:TruncationStrategy})
+    trunc_y = circshift(rotr90(trunc_x), (1, 0))
+    for (i, t) in enumerate(trunc_y)
+        if t isa TruncationSpace
+            trunc_y[i] = truncspace(flip(t.space)')
+        end
+    end
+    return trunc_y
+end
+function _rot180_trunc_x(trunc_x::AbstractMatrix{<:TruncationStrategy})
+    trunc_x_ = circshift(rot180(trunc_x), (0, -1))
+    for (i, t) in enumerate(trunc_x_)
+        trunc_x_[i] = truncspace(flip(t.space)')
+    end
+    return trunc_x_
+end
 
-    if directions == 2
-        truncs_rotated[NORTH, :, :] = circshift(
-            rotl90(trunc.truncs[EAST, :, :]), (0, -1)
-        )
-        truncs_rotated[EAST, :, :] = rotl90(trunc.truncs[NORTH, :, :])
-    elseif directions == 4
-        for dir in 1:4
-            dir′ = _prev(dir, 4)
-            truncs_rotated[dir′, :, :] = rotl90(trunc.truncs[dir, :, :])
+function _rotl90_trunc_y(trunc_y::AbstractMatrix{<:TruncationStrategy})
+    trunc_x = circshift(rotl90(trunc_y), (0, -1))
+    for (i, t) in enumerate(trunc_x)
+        if t isa TruncationSpace
+            trunc_x[i] = truncspace(flip(t.space)')
         end
-    else
-        throw(ArgumentError("Unsupported number of directions for rotl90: $directions"))
     end
-    return SiteDependentTruncation(truncs_rotated)
+    return trunc_x
+end
+function _rotr90_trunc_y(trunc_y::AbstractMatrix{<:TruncationStrategy})
+    trunc_x = rotr90(trunc_y)
+    return trunc_x
+end
+function _rot180_trunc_y(trunc_y::AbstractMatrix{<:TruncationStrategy})
+    trunc_y_ = circshift(rot180(trunc_y), (1, 0))
+    for (i, t) in enumerate(trunc_y_)
+        trunc_y_[i] = truncspace(flip(t.space)')
+    end
+    return trunc_y_
+end
+
+function Base.rotl90(trunc::SiteDependentTruncation)
+    trunc_y = _rotl90_trunc_x(trunc[1, :, :])
+    trunc_x = _rotl90_trunc_y(trunc[2, :, :])
+    trunc = stack((trunc_x, trunc_y); dims = 1)
+    return SiteDependentTruncation(trunc)
+end
+function Base.rotr90(trunc::SiteDependentTruncation)
+    trunc_y = _rotr90_trunc_x(trunc[1, :, :])
+    trunc_x = _rotr90_trunc_y(trunc[2, :, :])
+    trunc = stack((trunc_x, trunc_y); dims = 1)
+    return SiteDependentTruncation(trunc)
+end
+function Base.rot180(trunc::SiteDependentTruncation)
+    trunc_x = _rot180_trunc_x(trunc[1, :, :])
+    trunc_y = _rot180_trunc_y(trunc[2, :, :])
+    trunc = stack((trunc_x, trunc_y); dims = 1)
+    return SiteDependentTruncation(trunc)
 end
diff --git a/src/environments/bp_environments.jl b/src/environments/bp_environments.jl
index 765d31c0c..ea7ae9301 100644
--- a/src/environments/bp_environments.jl
+++ b/src/environments/bp_environments.jl
@@ -158,6 +158,15 @@ function eachcoordinate(x::BPEnv, dirs)
     return collect(Iterators.product(dirs, axes(x, 2), axes(x, 3)))
 end
 
+## Bipartite check
+function _is_bipartite(env::BPEnv)
+    (size(env, 2) == size(env, 3) == 2) || (return false)
+    for (d, c) in Iterators.product(axes(env, 1), axes(env, 3))
+        (env[d, 1, c] == env[d, 2, _next(c, 2)]) || (return false)
+    end
+    return true
+end
+
 # conversion to CTMRGEnv
 """
     CTMRGEnv(bp_env::BPEnv)
diff --git a/src/environments/suweight.jl b/src/environments/suweight.jl
index 2966d4155..cd24cb111 100644
--- a/src/environments/suweight.jl
+++ b/src/environments/suweight.jl
@@ -138,6 +138,15 @@ TensorKit.spacetype(::Type{T}) where {E, T <: SUWeight{E}} = spacetype(E)
 TensorKit.sectortype(w::SUWeight) = sectortype(typeof(w))
 TensorKit.sectortype(::Type{<:SUWeight{T}}) where {T} = sectortype(spacetype(T))
 
+## Bipartite check
+function _is_bipartite(wts::SUWeight)
+    (size(wts, 2) == size(wts, 3) == 2) || (return false)
+    for (d, c) in Iterators.product(1:2, 1:2)
+        (wts[d, 1, c] == wts[d, 2, _next(c, 2)]) || (return false)
+    end
+    return true
+end
+
 ## (Approximate) equality
 function Base.:(==)(wts1::SUWeight, wts2::SUWeight)
     return wts1.data == wts2.data
diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl
index 172c56b05..0eb44444a 100644
--- a/src/operators/infinitepepo.jl
+++ b/src/operators/infinitepepo.jl
@@ -202,6 +202,14 @@ function InfinitePEPS(ρ::InfinitePEPO)
     )
 end
 
+## Bipartite check for PEPS/PEPO
+function _is_bipartite(psi::InfiniteState)
+    (size(psi, 1) == size(psi, 2) == 2) || (return false)
+    for (c, h) in Iterators.product(1:2, 1:size(psi, 3))
+        (psi[1, c, h] == psi[2, _next(c, 2), h]) || (return false)
+    end
+    return true
+end
 
 ## Vector interface
 
diff --git a/test/bondenv/benv_ctm.jl b/test/bondenv/benv_ctm.jl
index 7676c5139..95ece1658 100644
--- a/test/bondenv/benv_ctm.jl
+++ b/test/bondenv/benv_ctm.jl
@@ -19,7 +19,7 @@ function get_hubbard_peps(t::Float64 = 1.0, U::Float64 = 8.0)
     wts = SUWeight(peps)
     alg = SimpleUpdate(; trunc = truncerror(; atol = 1.0e-10) & truncrank(4))
     evolver = TimeEvolver(peps, H, 1.0e-2, 10000, alg, wts)
-    peps, = time_evolve(evolver; tol = 1.0e-8, check_interval = 2000)
+    peps, = time_evolve(evolver, H; tol = 1.0e-8, check_interval = 2000)
     normalize!.(peps.A, Inf)
     return peps
 end
diff --git a/test/bp/gaugefix.jl b/test/bp/gaugefix.jl
index 941c8cbea..76d8fdc4e 100644
--- a/test/bp/gaugefix.jl
+++ b/test/bp/gaugefix.jl
@@ -3,54 +3,80 @@ using Random
 using TensorKit
 using PEPSKit
 using PEPSKit: compare_weights, random_dual!, twistdual
+using PEPSKit: _next, _is_bipartite
 
-@testset "Compare BP and SU ($S, posdef msgs = $h)" for (S, h) in
-    Iterators.product([U1Irrep, FermionParity], [true, false])
-    unitcell = (2, 3)
+@testset "BP vs SU ($S, bipartite = $(bipartite), posdef msgs = $h)" for
+    (S, bipartite, h) in Iterators.product(
+        [U1Irrep, FermionParity], [true, false], [true, false]
+    )
+    unitcell = bipartite ? (2, 2) : (2, 3)
     elt = ComplexF64
     maxiter, tol = 100, 1.0e-9
     Random.seed!(52840679)
     Pspaces, Nspaces, Espaces = if S == U1Irrep
-        map(zip(rand(1:2, unitcell), rand(1:2, unitcell), rand(1:2, unitcell))) do (d0, d1, d2)
+        map(rand(1:2, unitcell), rand(1:2, unitcell), rand(1:2, unitcell)) do d0, d1, d2
                 Vect[S](0 => d0, 1 => d1, -1 => d2)
         end,
-            map(zip(rand(2:4, unitcell), rand(2:4, unitcell), rand(2:4, unitcell))) do (d0, d1, d2)
+            map(rand(2:4, unitcell), rand(2:4, unitcell), rand(2:4, unitcell)) do d0, d1, d2
                 Vect[S](0 => d0, 1 => d1, -1 => d2)
         end,
-            map(zip(rand(2:4, unitcell), rand(2:4, unitcell), rand(2:4, unitcell))) do (d0, d1, d2)
+            map(rand(2:4, unitcell), rand(2:4, unitcell), rand(2:4, unitcell)) do d0, d1, d2
                 Vect[S](0 => d0, 1 => d1, -1 => d2)
         end
     else
-        map(zip(rand(2:3, unitcell), rand(2:3, unitcell))) do (d0, d1)
+        map(rand(2:3, unitcell), rand(2:3, unitcell)) do d0, d1
                 Vect[S](0 => d0, 1 => d1)
         end,
-            map(zip(rand(2:4, unitcell), rand(2:4, unitcell))) do (d0, d1)
+            map(rand(2:4, unitcell), rand(2:4, unitcell)) do d0, d1
                 Vect[S](0 => d0, 1 => d1)
         end,
-            map(zip(rand(2:4, unitcell), rand(2:4, unitcell))) do (d0, d1)
+            map(rand(2:4, unitcell), rand(2:4, unitcell)) do d0, d1
                 Vect[S](0 => d0, 1 => d1)
         end
     end
     Nspaces, Espaces = random_dual!(Nspaces), random_dual!(Espaces)
+    if bipartite
+        for c in 1:2
+            cp1 = _next(c, 2)
+            Pspaces[2, c] = Pspaces[1, cp1]
+            Nspaces[2, c] = Nspaces[1, cp1]
+            Espaces[2, c] = Espaces[1, cp1]
+        end
+    end
     peps0 = InfinitePEPS(randn, elt, Pspaces, Nspaces, Espaces)
+    if bipartite
+        for c in 1:2
+            peps0[2, c] = copy(peps0[1, _next(c, 2)])
+        end
+    end
 
     # start by gauging with SU
     peps1, wts1 = gauge_fix(peps0, SUGauge(; maxiter, tol))
     for (a0, a1) in zip(peps0.A, peps1.A)
         @test space(a0) == space(a1)
     end
+    if bipartite
+        @test _is_bipartite(peps1)
+        @test _is_bipartite(wts1)
+    end
     normalize!.(wts1.data)
 
     # find BP fixed point and SUWeight
-    bp_alg = BeliefPropagation(; maxiter, tol, project_hermitian = h)
+    bp_alg = BeliefPropagation(; maxiter, tol, bipartite, project_hermitian = h)
     env = BPEnv(randn, elt, peps1; posdef = h)
     env, err = leading_boundary(env, peps1, bp_alg)
+    if bipartite
+        @test _is_bipartite(env)
+    end
     wts2 = SUWeight(env)
     normalize!.(wts2.data)
     @test compare_weights(wts1, wts2) < 1.0e-9
 
     bpg_alg = BPGauge()
     peps2, XXinv = @constinferred gauge_fix(peps1, bpg_alg, env)
+    if bipartite
+        @test _is_bipartite(peps2)
+    end
     for (a1, a2) in zip(peps1.A, peps2.A)
         @test space(a1) == space(a2)
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index 69f670e35..a3f760cff 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -104,7 +104,7 @@ end
         @time @safetestset "Cluster truncation with projectors" begin
             include("timeevol/cluster_projectors.jl")
         end
-        @time @safetestset "Time evolution with site-dependent truncation" begin
+        @time @safetestset "Fixed-space and site-dependent truncation" begin
             include("timeevol/sitedep_truncation.jl")
         end
         @time @safetestset "Transverse field Ising model at finite temperature" begin
diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl
index 53fad0a68..1276ad2f8 100644
--- a/test/timeevol/cluster_projectors.jl
+++ b/test/timeevol/cluster_projectors.jl
@@ -35,7 +35,8 @@ Vspaces = [
         flips = [isdual(space(M, 1)) for M in Ms1[2:end]]
         # no truncation
         Ms2 = _flip_virtuals!(deepcopy(Ms1), flips)
-        wts2, ϵs, = _cluster_truncate!(Ms2, fill(FixedSpaceTruncation(), N - 1))
+        truncs = [truncrank(dim(space(M, 1))) for M in Ms2[2:end]]
+        wts2, ϵs, = _cluster_truncate!(Ms2, truncs)
         @test all((ϵ == 0) for ϵ in ϵs)
         normalize!.(Ms2, Inf)
         @test fidelity_cluster(Ms1, Ms2) ≈ 1.0
diff --git a/test/timeevol/sitedep_truncation.jl b/test/timeevol/sitedep_truncation.jl
index a36f3de75..36952f001 100644
--- a/test/timeevol/sitedep_truncation.jl
+++ b/test/timeevol/sitedep_truncation.jl
@@ -1,70 +1,63 @@
 using Test
-using LinearAlgebra
 using Random
 using TensorKit
 using PEPSKit
-using PEPSKit: NORTH, EAST, _next
+using PEPSKit: _is_bipartite, _get_fixedspacetrunc
 
-function get_bonddims(peps::InfinitePEPS)
-    xdims = collect(dim(domain(t, EAST)) for t in peps.A)
-    ydims = collect(dim(domain(t, NORTH)) for t in peps.A)
-    return stack([xdims, ydims]; dims = 1)
-end
-
-function get_bonddims(wts::SUWeight)
-    xdims = collect(dim(space(wt, 1)) for wt in wts[1, :, :])
-    ydims = collect(dim(space(wt, 1)) for wt in wts[2, :, :])
-    return stack([xdims, ydims]; dims = 1)
-end
+elt = Float64
+Nr, Nc = 2, 2
+Vps = fill(U1Space(1 / 2 => 1, -1 / 2 => 1), (Nr, Nc))
+Vns = [
+    U1Space(0 => 1, 1 => 2, -1 => 1) U1Space(0 => 1, 1 => 2, -1 => 1)';
+    U1Space(0 => 1, 1 => 2, -1 => 1)' U1Space(0 => 1, 1 => 2, -1 => 1)
+]
+Ves1 = [
+    U1Space(1 / 2 => 1, -1 / 2 => 2, -3 / 2 => 1)' U1Space(0 => 1, 1 => 1, -1 => 2);
+    U1Space(0 => 1, 1 => 1, -1 => 2) U1Space(1 / 2 => 1, -1 / 2 => 2, -3 / 2 => 1)'
+]
+Ves2 = [
+    U1Space(0 => 1, 1 => 2, -1 => 1)' U1Space(0 => 1, 1 => 1, -1 => 2);
+    U1Space(0 => 1, 1 => 1, -1 => 2) U1Space(0 => 1, 1 => 2, -1 => 1)'
+]
+Venv = U1Space(0 => 2, 1 => 1, -1 => 1)
+Random.seed!(48736)
+states = (
+    InfinitePEPS(randn, elt, Vps, Vns, Ves1),
+    InfinitePEPO(randn, elt, Vps, Vns, Ves2),
+)
 
-@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))
-    Random.seed!(100)
-    peps0 = InfinitePEPS(rand, Float64, ℂ^2, ℂ^10; unitcell = (Nr, Nc))
-    # make state bipartite
-    for r in 1:2
-        peps0.A[_next(r, 2), 2] = copy(peps0.A[r, 1])
-    end
-    env0 = SUWeight(peps0)
-    normalize!.(peps0.A, Inf)
-    # 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(; 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
-    for col in 1:2
-        cp1 = PEPSKit._next(col, 2)
-        @test (
-            peps.A[1, col] == peps.A[2, cp1] &&
-                env[1, 1, col] == env[1, 2, cp1] &&
-                env[2, 1, col] == env[2, 2, cp1]
+@testset "Rotation of SiteDependentTruncation" begin
+    state = states[1]
+    for f in (rotl90, rotr90, rot180)
+        trunc1 = f(_get_fixedspacetrunc(state))
+        trunc2 = _get_fixedspacetrunc(f(state))
+        @test all(
+            t1.space == t2.space for (t1, t2) in zip(trunc1.truncs, trunc2.truncs)
         )
     end
 end
 
-@testset "Simple update: generic 2-site and 3-site" begin
-    Nr, Nc = 3, 4
-    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
-    trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims))
-    alg = SimpleUpdate(; trunc)
-    # 2-site SU
-    ham = j1_j2_model(Float64, Trivial, InfiniteSquare(Nr, Nc); J2 = 0.0, sublattice = false)
-    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
-    ham = j1_j2_model(Float64, Trivial, InfiniteSquare(Nr, Nc); J2 = 0.2, sublattice = false)
-    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
-    @test get_bonddims(peps) == bonddims
-    @test get_bonddims(env) == bonddims
+@testset "Simple update on $(typeof(state0).name.wrapper), bipartite = $(bipartite)" for
+    (state0, bipartite) in Iterators.product(states, (true, false))
+    J2 = 0.5
+    if bipartite
+        state0[2, 1] = copy(state0[1, 2])
+        state0[2, 2] = copy(state0[1, 1])
+        J2 = 0.0
+    end
+    ham = j1_j2_model(elt, U1Irrep, InfiniteSquare(Nr, Nc); J1 = 1.0, J2, sublattice = false)
+    # converted internally to SiteDependentTruncation
+    alg = SimpleUpdate(; trunc = FixedSpaceTruncation(), bipartite)
+    wts0 = SUWeight(state0)
+    state, wts, = time_evolve(state0, ham, 0.1, 1, alg, wts0)
+    for (t, t0) in zip(state.A, state0.A)
+        @test space(t) == space(t0)
+    end
+    for (wt, wt0) in zip(wts.data, wts0.data)
+        @test space(wt) == space(wt0)
+    end
+    if bipartite
+        @test _is_bipartite(state)
+        @test _is_bipartite(wts)
+    end
 end