diff --git a/src/algorithms/time_evolution/evoltools.jl b/src/algorithms/time_evolution/evoltools.jl
index 588fd9c14..4dfab91da 100644
--- a/src/algorithms/time_evolution/evoltools.jl
+++ b/src/algorithms/time_evolution/evoltools.jl
@@ -1,11 +1,11 @@
 """
-    get_expham(dt::Number, H::LocalOperator)
+    get_expham(H::LocalOperator, dt::Number)
 
 Compute `exp(-dt * op)` for each term `op` in `H`,
 and combine them into a new LocalOperator.
 Each `op` in `H` must be a single `TensorMap`.
 """
-function get_expham(dt::Number, H::LocalOperator)
+function get_expham(H::LocalOperator, dt::Number)
     return LocalOperator(
         physicalspace(H), (sites => exp(-dt * op) for (sites, op) in H.terms)...
     )
@@ -181,19 +181,30 @@ Obtain the 3-site gate MPO on the southeast cluster at position `[row, col]`
         c      c+1
 ```
 """
-function _get_gatempo_se(gate::LocalOperator, row::Int, col::Int)
-    Nr, Nc = size(gate.lattice)
+function _get_gatempo_se(ham::LocalOperator, dt::Number, row::Int, col::Int)
+    Nr, Nc = size(ham)
     @assert 1 <= row <= Nr && 1 <= col <= Nc
-    unit = id(space(gate.terms[1].second, 1))
-    sites = (
+    sites = [
         CartesianIndex(row, col),
         CartesianIndex(row, col + 1),
         CartesianIndex(row - 1, col + 1),
-    )
-    nb1x = get_gateterm(gate, (sites[1], sites[2]))
-    nb1y = get_gateterm(gate, (sites[2], sites[3]))
-    nb2 = get_gateterm(gate, (sites[1], sites[3]))
-    op = (1 / 2) * (nb1x ⊗ unit + unit ⊗ nb1y) + permute(nb2 ⊗ unit, ((1, 3, 2), (4, 6, 5)))
+    ]
+    nb1x = get_gateterm(ham, (sites[1], sites[2]))
+    nb1y = get_gateterm(ham, (sites[2], sites[3]))
+    nb2 = get_gateterm(ham, (sites[1], sites[3]))
+    # identity operator at each site
+    units = map(sites) do site
+        site_ = CartesianIndex(mod1(site[1], Nr), mod1(site[2], Nc))
+        return id(physicalspace(ham)[site_])
+    end
+    # when iterating through ┘, └, ┌, ┐ clusters in the unit cell,
+    # NN / NNN bonds are counted 4 / 2 times, respectively.
+    @tensor Odt[i' j' k'; i j k] :=
+        -dt * (
+            (nb1x[i' j'; i j] * units[3][k' k] + units[1][i'; i] * nb1y[j' k'; j k]) / 4 +
+            (nb2[i' k'; i k] * units[2][j'; j]) / 2
+        )
+    op = exp(Odt)
     return gate_to_mpo3(op)
 end
 
@@ -201,7 +212,7 @@ end
 Construct the 3-site gate MPOs on the southeast cluster 
 for 3-site simple update on square lattice.
 """
-function _get_gatempos_se(gate::LocalOperator)
-    Nr, Nc = size(gate.lattice)
-    return collect(_get_gatempo_se(gate, r, c) for r in 1:Nr, c in 1:Nc)
+function _get_gatempos_se(ham::LocalOperator, dt::Number)
+    Nr, Nc = size(ham.lattice)
+    return collect(_get_gatempo_se(ham, dt, r, c) for r in 1:Nr, c in 1:Nc)
 end
diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index de75c323d..01eb070d0 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -9,7 +9,7 @@ Each SU run is converged when the singular value difference becomes smaller than
 $(TYPEDFIELDS)
 """
 struct SimpleUpdate
-    dt::Float64
+    dt::Number
     tol::Float64
     maxiter::Int
     trscheme::TensorKit.TruncationScheme
@@ -29,7 +29,7 @@ Simple update of the x-bond `peps.weights[1,r,c]`.
                 [2,r+1,c]           [2,r+1,c+1]
 ```
 """
-function _su_bondx!(
+function _su_xbond!(
     row::Int,
     col::Int,
     gate::AbstractTensorMap{T,S,2,2},
@@ -94,7 +94,7 @@ function su_iter(
                     direction == 1 ? gate : gate_mirrored,
                     (CartesianIndex(r, 1), CartesianIndex(r, 2)),
                 )
-                ϵ = _su_bondx!(r, 1, term, peps2, alg)
+                ϵ = _su_xbond!(r, 1, term, peps2, alg)
                 peps2.vertices[rp1, 2] = deepcopy(peps2.vertices[r, 1])
                 peps2.vertices[rp1, 1] = deepcopy(peps2.vertices[r, 2])
                 peps2.weights[1, rp1, 2] = deepcopy(peps2.weights[1, r, 1])
@@ -106,7 +106,7 @@ function su_iter(
                     direction == 1 ? gate : gate_mirrored,
                     (CartesianIndex(r, c), CartesianIndex(r, c + 1)),
                 )
-                ϵ = _su_bondx!(r, c, term, peps2, alg)
+                ϵ = _su_xbond!(r, c, term, peps2, alg)
             end
         end
         if direction == 2
@@ -128,7 +128,7 @@ function _simpleupdate2site(
 )
     time_start = time()
     # exponentiating the 2-site Hamiltonian gate
-    gate = get_expham(alg.dt, ham)
+    gate = get_expham(ham, alg.dt)
     wtdiff = 1.0
     wts0 = deepcopy(peps.weights)
     for count in 1:(alg.maxiter)
diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl
index 678ac69a0..1ce887b5b 100644
--- a/src/algorithms/time_evolution/simpleupdate3site.jl
+++ b/src/algorithms/time_evolution/simpleupdate3site.jl
@@ -414,13 +414,13 @@ function su3site_iter(
         ),
     )
     peps2 = deepcopy(peps)
-    for i in 1:2
+    for i in 1:4
         for site in CartesianIndices(peps2.vertices)
             r, c = site[1], site[2]
             gs = gatempos[i][r, c]
             _su3site_se!(r, c, gs, peps2, alg)
         end
-        peps2 = (i == 1) ? rotl90(peps2) : rotr90(peps2)
+        peps2 = rotl90(peps2)
     end
     return peps2
 end
@@ -432,9 +432,13 @@ function _simpleupdate3site(
     peps::InfiniteWeightPEPS, ham::LocalOperator, alg::SimpleUpdate; check_interval::Int=500
 )
     time_start = time()
-    gate = get_expham(alg.dt, ham)
-    # convert gates to 3-site MPOs
-    gatempos = [_get_gatempos_se(gate), _get_gatempos_se(rotl90(gate))]
+    # convert Hamiltonian to 3-site exponentiated gate MPOs
+    gatempos = [
+        _get_gatempos_se(ham, alg.dt),
+        _get_gatempos_se(rotl90(ham), alg.dt),
+        _get_gatempos_se(rot180(ham), alg.dt),
+        _get_gatempos_se(rotr90(ham), alg.dt),
+    ]
     wtdiff = 1.0
     wts0 = deepcopy(peps.weights)
     for count in 1:(alg.maxiter)
diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl
index 0bafec669..1d925792e 100644
--- a/src/operators/localoperator.jl
+++ b/src/operators/localoperator.jl
@@ -105,6 +105,8 @@ function physicalspace(O::LocalOperator)
     return O.lattice
 end
 
+Base.size(O::LocalOperator) = size(physicalspace(O))
+
 # Real and imaginary part
 # -----------------------
 function Base.real(O::LocalOperator)