Fix 3-site SU dt and reduce artificial C4v breaking

src/algorithms/time_evolution/evoltools.jl

@@ -181,27 +181,38 @@
         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
 
 """
 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

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

src/algorithms/time_evolution/simpleupdate3site.jl

@@ -414,13 +414,13 @@
         ),
     )
     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 @@
     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)

src/operators/localoperator.jl

@@ -105,6 +105,10 @@
     return O.lattice
 end
 
+function Base.size(O::LocalOperator)
+    return size(O.lattice)
+end
+
 # Real and imaginary part
 # -----------------------
 function Base.real(O::LocalOperator)