diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl
index 2377561a..b7b00467 100644
--- a/src/algorithms/ctmrg/ctmrg.jl
+++ b/src/algorithms/ctmrg/ctmrg.jl
@@ -137,7 +137,7 @@ end
 Perform CTMRG left move on the `col`-th column
 """
 function ctmrg_leftmove(col::Int, state, envs::CTMRGEnv, alg::SequentialCTMRG)
-    """
+    #=
         ----> left move
         C1 ← T1 ←   r-1
         ↓    ‖
@@ -145,7 +145,7 @@ function ctmrg_leftmove(col::Int, state, envs::CTMRGEnv, alg::SequentialCTMRG)
         ↓    ‖
         C4 → T3 →   r+1
         c-1  c 
-    """
+    =#
     projectors, info = ctmrg_projectors(col, state, envs, alg)
     envs = ctmrg_renormalize(col, projectors, state, envs, alg)
     return envs, info
diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl
index a69f9bb4..64b80984 100644
--- a/src/operators/localoperator.jl
+++ b/src/operators/localoperator.jl
@@ -187,4 +187,3 @@ function Base.rot180(H::LocalOperator)
     )
     return LocalOperator(lattice2, terms2...)
 end
-
diff --git a/src/utility/util.jl b/src/utility/util.jl
index 68ef13c0..66218978 100644
--- a/src/utility/util.jl
+++ b/src/utility/util.jl
@@ -41,11 +41,12 @@ function ChainRulesCore.rrule(
     ::typeof(sdiag_pow), S::AbstractTensorMap, pow::Real; tol::Real=eps(eltype(S))^(3 / 4)
 )
     tol *= norm(S, Inf)
-    spow = sdiag_pow(S, pow - 1; tol)
+    spow = sdiag_pow(S, pow; tol)
+    spow2 = sdiag_pow(S, pow - 1; tol)
     function sdiag_pow_pullback(c̄)
-        return (ChainRulesCore.NoTangent(), pow * _elementwise_mult(c̄, spow))
+        return (ChainRulesCore.NoTangent(), pow * _elementwise_mult(c̄, spow2))
     end
-    return invsq, sdiag_pow_pullback
+    return spow, sdiag_pow_pullback
 end
 
 # Check whether diagonals contain degenerate values up to absolute or relative tolerance
diff --git a/test/heisenberg_suad.jl b/test/heisenberg_suad.jl
deleted file mode 100644
index 7902c989..00000000
--- a/test/heisenberg_suad.jl
+++ /dev/null
@@ -1,78 +0,0 @@
-using Test
-using Printf
-using Random
-using PEPSKit
-using TensorKit
-using OptimKit
-using KrylovKit
-import Statistics: mean
-import MPSKitModels: S_x, S_y, S_z, S_exchange
-include("utility/measure_heis.jl")
-import .MeasureHeis: measure_heis
-
-# random initialization of 2x2 iPEPS with weights and CTMRGEnv (using real numbers)
-Dcut, χenv = 4, 16
-N1, N2 = 2, 2
-Random.seed!(0)
-peps = InfiniteWeightPEPS(rand, Float64, ℂ^2, ℂ^Dcut; unitcell=(N1, N2))
-# normalize vertex tensors
-for ind in CartesianIndices(peps.vertices)
-    peps.vertices[ind] /= norm(peps.vertices[ind], Inf)
-end
-
-# Heisenberg model Hamiltonian
-# (only includes nearest neighbor terms)
-lattice = InfiniteSquare(N1, N2)
-onsite = TensorMap([1.0 0.0; 0.0 1.0], ℂ^2, ℂ^2)
-ham = heisenberg_XYZ(lattice; Jx=1.0, Jy=1.0, Jz=1.0)
-# convert to real tensors
-ham = LocalOperator(ham.lattice, Tuple(ind => real(op) for (ind, op) in ham.terms)...)
-
-# Include the onsite operators in two ways
-ham_SU = LocalOperator(
-    ham.lattice,
-    Tuple(
-        sites => op + (S_z() ⊗ onsite) / 2 for
-        (sites, op) in ham.terms if length(sites) == 2
-    )...,
-)
-ham_CTMRG = LocalOperator(
-    ham.lattice,
-    Tuple(ind => op for (ind, op) in ham.terms)...,
-    ((idx,) => S_z() for idx in vertices(lattice))...,
-)
-
-# simple update with ham_SU
-dts = [1e-2, 1e-3, 4e-4, 1e-4]
-tols = [1e-6, 1e-8, 1e-8, 1e-8]
-maxiter = 10000
-for (n, (dt, tol)) in enumerate(zip(dts, tols))
-    Dcut2 = (n == 1 ? Dcut + 1 : Dcut)
-    trscheme = truncerr(1e-10) & truncdim(Dcut2)
-    alg = SimpleUpdate(dt, tol, maxiter, trscheme)
-    result = simpleupdate(peps, ham_SU, alg; bipartite=false)
-    global peps = result[1]
-end
-# absort weight into site tensors
-peps = InfinitePEPS(peps)
-# CTMRG
-envs = CTMRGEnv(rand, Float64, peps, ℂ^χenv)
-trscheme = truncerr(1e-9) & truncdim(χenv)
-ctm_alg = CTMRG(; tol=1e-10, verbosity=2, trscheme=trscheme, ctmrgscheme=:simultaneous)
-envs = leading_boundary(envs, peps, ctm_alg)
-# measure physical quantities
-meas = measure_heis(peps, ham_SU, envs)
-display(meas)
-
-# continue with auto-diff optimization
-opt_alg = PEPSOptimize(;
-    boundary_alg=ctm_alg,
-    optimizer=LBFGS(4; maxiter=100, gradtol=1e-3, verbosity=2),
-    gradient_alg=LinSolver(; solver=GMRES(; tol=1e-6), iterscheme=:diffgauge),
-    reuse_env=true,
-)
-result = fixedpoint(peps, ham_CTMRG, opt_alg, envs)
-meas2 = measure_heis(result.peps, ham_CTMRG, result.env)
-display(meas2)
-
-@test isapprox(result.E / (N1 * N2), meas["e_site"], atol=1e-2)
diff --git a/test/heisenberg_sufu.jl b/test/heisenberg_sufu.jl
index 36fc9706..32b3095e 100644
--- a/test/heisenberg_sufu.jl
+++ b/test/heisenberg_sufu.jl
@@ -4,8 +4,6 @@ using Random
 using PEPSKit
 using TensorKit
 import Statistics: mean
-include("utility/measure_heis.jl")
-import .MeasureHeis: measure_heis
 
 # benchmark data is from Phys. Rev. B 94, 035133 (2016)
 
@@ -39,7 +37,7 @@ ham = LocalOperator(ham.lattice, Tuple(ind => real(op) for (ind, op) in ham.term
 # simple update
 dts = [1e-2, 1e-3, 4e-4, 1e-4]
 tols = [1e-6, 1e-8, 1e-8, 1e-8]
-maxiter = 10000
+maxiter = 5000
 for (n, (dt, tol)) in enumerate(zip(dts, tols))
     trscheme = truncerr(1e-10) & truncdim(Dcut)
     alg = SimpleUpdate(dt, tol, maxiter, trscheme)
@@ -54,9 +52,6 @@ trscheme = truncerr(1e-9) & truncdim(χenv)
 ctm_alg = CTMRG(; tol=1e-10, verbosity=2, trscheme=trscheme, ctmrgscheme=:sequential)
 envs = leading_boundary(envs, peps, ctm_alg)
 # measure physical quantities
-meas = measure_heis(peps, ham, envs)
-display(meas)
-@info @sprintf("Energy = %.8f\n", meas["e_site"])
-@info @sprintf("Staggered magnetization = %.8f\n", mean(meas["mag_norm"]))
-@test isapprox(meas["e_site"], -0.6675; atol=1e-3)
-@test isapprox(mean(meas["mag_norm"]), 0.3767; atol=1e-3)
+e_site = costfun(peps, envs, ham)
+@info @sprintf("Energy = %.8f\n", e_site)
+@test isapprox(e_site, -0.6675; atol=1e-3)
diff --git a/test/runtests.jl b/test/runtests.jl
index 45bd68e3..19008ebb 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -56,9 +56,6 @@ end
         @time @safetestset "Heisenberg model (simple and full update)" begin
             include("heisenberg_sufu.jl")
         end
-        @time @safetestset "Hubbard model (simple update)" begin
-            include("hubbard_su.jl")
-        end
         @time @safetestset "J1-J2 model" begin
             include("j1j2_model.jl")
         end
diff --git a/test/utility/measure_heis.jl b/test/utility/measure_heis.jl
deleted file mode 100644
index e63b1b86..00000000
--- a/test/utility/measure_heis.jl
+++ /dev/null
@@ -1,50 +0,0 @@
-module MeasureHeis
-
-export measure_heis
-
-using TensorKit
-import MPSKitModels: S_x, S_y, S_z, S_exchange
-using PEPSKit
-using Statistics: mean
-
-"""
-Measure magnetization on each site
-"""
-function cal_mags(peps::InfinitePEPS, envs::CTMRGEnv)
-    Nr, Nc = size(peps)
-    lattice = collect(space(t, 1) for t in peps.A)
-    # detect symmetry on physical axis
-    symm = sectortype(space(peps.A[1, 1]))
-    if symm == Trivial
-        Sas = real.([S_x(symm), im * S_y(symm), S_z(symm)])
-    elseif symm == U1Irrep
-        # only Sz preserves <Sz>
-        Sas = real.([S_z(symm)])
-    else
-        throw(ArgumentError("Unrecognized symmetry on physical axis"))
-    end
-    return [
-        collect(
-            expectation_value(
-                peps, LocalOperator(lattice, (CartesianIndex(r, c),) => Sa), envs
-            ) for (r, c) in Iterators.product(1:Nr, 1:Nc)
-        ) for Sa in Sas
-    ]
-end
-
-"""
-Measure physical quantities for Heisenberg model
-"""
-function measure_heis(peps::InfinitePEPS, H::LocalOperator, envs::CTMRGEnv)
-    results = Dict{String,Any}()
-    Nr, Nc = size(peps)
-    results["e_site"] = costfun(peps, envs, H) / (Nr * Nc)
-    results["mag"] = cal_mags(peps, envs)
-    results["mag_norm"] = collect(
-        norm([mags[r, c] for mags in results["mag"]]) for
-        (r, c) in Iterators.product(1:Nr, 1:Nc)
-    )
-    return results
-end
-
-end