diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml
index c88bccd70..e03fc30b0 100644
--- a/.github/workflows/Tests.yml
+++ b/.github/workflows/Tests.yml
@@ -1,51 +1,26 @@
-name: Tests
+name: CI
+
 on:
   push:
     branches:
-      - 'master'
       - 'main'
-      - 'release-'
     tags: '*'
     paths-ignore:
       - 'docs/**'
   pull_request:
+    types: [opened, synchronize, reopened, ready_for_review, converted_to_draft]
   workflow_dispatch:
 
 concurrency:
   group: ${{ github.workflow }}-${{ github.ref }}
-  # Cancel intermediate builds: only if it is a pull request build.
   cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
 
 jobs:
   tests:
     name: "Tests"
-    strategy:
-      fail-fast: false
-      matrix:
-        version:
-          - 'lts' # minimal supported version
-          - '1'
-        group:
-          - bp
-          - types
-          - ctmrg
-          - boundarymps
-          - gradients
-          - examples
-          - utility
-          - bondenv
-          - timeevol
-          - toolbox
-        os:
-          - ubuntu-latest
-          - macOS-latest
-          - windows-latest
-    uses: "QuantumKitHub/QuantumKitHubActions/.github/workflows/Tests.yml@main"
+    uses: "QuantumKitHub/QuantumKitHubActions/.github/workflows/TestGroups.yml@main"
     with:
-      group: "${{ matrix.group }}"
-      julia-version: "${{ matrix.version }}"
-      os: "${{ matrix.os }}"
+      fast: "${{ github.event.pull_request.draft == true }}"
       timeout-minutes: 120
     secrets:
       CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
-
diff --git a/Project.toml b/Project.toml
index b345da7af..6e16222f8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "PEPSKit"
 uuid = "52969e89-939e-4361-9b68-9bc7cde4bdeb"
-authors = ["Paul Brehmer", "Lander Burgelman", "Zhengyuan Yue", "Lukas Devos <ldevos98@gmail.com>"]
 version = "0.7.0"
+authors = ["Paul Brehmer", "Lander Burgelman", "Zhengyuan Yue", "Lukas Devos <ldevos98@gmail.com>"]
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
@@ -28,11 +28,8 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 [compat]
 Accessors = "0.1"
 ChainRulesCore = "1.0"
-ChainRulesTestUtils = "1.13"
-SafeTestsets = "0.1"
 Compat = "3.46, 4.2"
 DocStringExtensions = "0.9.3"
-FiniteDifferences = "0.12"
 KrylovKit = "0.9.5, 0.10"
 LinearAlgebra = "1"
 LoggingExtras = "1"
@@ -42,23 +39,13 @@ MatrixAlgebraKit = "0.6.5"
 OhMyThreads = "0.7, 0.8"
 OptimKit = "0.4"
 Printf = "1"
-QuadGK = "2.11.1"
 Random = "1"
 Statistics = "1"
 TensorKit = "0.16.2"
 TensorOperations = "5"
-TestExtras = "0.3"
 VectorInterface = "0.4, 0.5"
 Zygote = "0.6, 0.7"
 julia = "1.10"
 
-[extras]
-ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
-FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
-QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
-
-[targets]
-test = ["Test", "TestExtras", "SafeTestsets", "ChainRulesTestUtils", "FiniteDifferences", "QuadGK"]
+[workspace]
+projects = ["test", "docs"]
diff --git a/test/Project.toml b/test/Project.toml
new file mode 100644
index 000000000..a64eeae63
--- /dev/null
+++ b/test/Project.toml
@@ -0,0 +1,32 @@
+name = "PEPSKitTests"
+
+[deps]
+Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
+ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
+FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
+KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
+MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502"
+MPSKitModels = "ca635005-6f8c-4cd1-b51d-8491250ef2ab"
+OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
+ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
+PEPSKit = "52969e89-939e-4361-9b68-9bc7cde4bdeb"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
+
+[compat]
+ChainRulesTestUtils = "1.13"
+FiniteDifferences = "0.12"
+ParallelTestRunner = "2.6.0"
+QuadGK = "2.11.1"
+Test = "1"
+TestExtras = "0.3"
+
+[sources]
+PEPSKit = {path = ".."}
diff --git a/test/ctmrg/contractions.jl b/test/ctmrg/contractions.jl
index ea589b3ee..7a91de3c2 100644
--- a/test/ctmrg/contractions.jl
+++ b/test/ctmrg/contractions.jl
@@ -2,7 +2,6 @@ using Test
 using Random
 using PEPSKit
 using TensorKit
-using TensorOperations: tensorcontract
 
 using PEPSKit: eachcoordinate, _next_coordinate
 using PEPSKit: EnlargedCorner, HalfInfiniteEnv, FullInfiniteEnv
diff --git a/test/runtests.jl b/test/runtests.jl
index a3f760cff..2283d9c89 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,159 +1,12 @@
-using SafeTestsets
+using ParallelTestRunner
+using PEPSKit
 
-# check if user supplied args
-pat = r"(?:--group=)(\w+)"
-arg_id = findfirst(contains(pat), ARGS)
-const GROUP = if isnothing(arg_id)
-    uppercase(get(ENV, "GROUP", "ALL"))
-else
-    uppercase(only(match(pat, ARGS[arg_id]).captures))
-end
+# --fast to indicate a smaller set of tests
+args = parse_args(ARGS; custom = ["fast"])
+fast = !isnothing(args.custom["fast"])
 
-@time begin
-    if GROUP == "ALL" || GROUP == "TYPES"
-        @time @safetestset "InfiniteSquareNetwork" begin
-            include("types/infinitesquarenetwork.jl")
-        end
-        @time @safetestset "InfinitePartitionFunction" begin
-            include("types/infinitepartitionfunction.jl")
-        end
-        @time @safetestset "SUWeight" begin
-            include("types/suweight.jl")
-        end
-        @time @safetestset "LocalOperator" begin
-            include("types/localoperator.jl")
-        end
-        @time @safetestset "LocalCircuit" begin
-            include("types/localcircuit.jl")
-        end
-    end
-    if GROUP == "ALL" || GROUP == "CTMRG"
-        @time @safetestset "Gauge Fixing" begin
-            include("ctmrg/gaugefix.jl")
-        end
-        @time @safetestset "Unit cell" begin
-            include("ctmrg/unitcell.jl")
-        end
-        @time @safetestset ":fixed CTMRG iteration scheme" begin
-            include("ctmrg/fixed_iterscheme.jl")
-        end
-        @time @safetestset "SUWeight conversion" begin
-            include("ctmrg/suweight.jl")
-        end
-        @time @safetestset "Flavors" begin
-            include("ctmrg/flavors.jl")
-        end
-        @time @safetestset "Jacobian real linearity" begin
-            include("ctmrg/jacobian_real_linear.jl")
-        end
-        @time @safetestset "Partition function" begin
-            include("ctmrg/partition_function.jl")
-        end
-        @time @safetestset "PEPO" begin
-            include("ctmrg/pepo.jl")
-        end
-        @time @safetestset "correlation length" begin
-            include("ctmrg/correlation_length.jl")
-        end
-        @time @safetestset "Contractions" begin
-            include("ctmrg/contractions.jl")
-        end
-    end
-    if GROUP == "ALL" || GROUP == "BP"
-        @time @safetestset "Unit cell bond matching" begin
-            include("bp/unitcell.jl")
-        end
-        @time @safetestset "Expectation values" begin
-            include("bp/expvals.jl")
-        end
-        @time @safetestset "Rotation of BPEnv" begin
-            include("bp/rotation.jl")
-        end
-        @time @safetestset "Gauge-fixing iPEPS" begin
-            include("bp/gaugefix.jl")
-        end
-    end
-    if GROUP == "ALL" || GROUP == "GRADIENTS"
-        @time @safetestset "CTMRG gradients" begin
-            include("gradients/ctmrg_gradients.jl")
-        end
-        @time @safetestset "C4v CTMRG gradients" begin
-            include("gradients/c4v_ctmrg_gradients.jl")
-        end
-    end
-    if GROUP == "ALL" || GROUP == "BOUNDARYMPS"
-        @time @safetestset "VUMPS" begin
-            include("boundarymps/vumps.jl")
-        end
-    end
-    if GROUP == "ALL" || GROUP == "BONDENV"
-        @time @safetestset "Iterative optimization after truncation" begin
-            include("bondenv/bond_truncate.jl")
-        end
-        @time @safetestset "Gauge fixing" begin
-            include("bondenv/benv_gaugefix.jl")
-        end
-        @time @safetestset "Full update bond environment" begin
-            include("bondenv/benv_ctm.jl")
-        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
-        @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
-            include("timeevol/tf_ising_finiteT.jl")
-        end
-        @time @safetestset "J1-J2 model at finite temperature" begin
-            include("timeevol/j1j2_finiteT.jl")
-        end
-    end
-    if GROUP == "ALL" || GROUP == "TOOLBOX"
-        @time @safetestset "Density matrix from double-layer PEPO" begin
-            include("toolbox/densitymatrices.jl")
-        end
-    end
-    if GROUP == "ALL" || GROUP == "UTILITY"
-        @time @safetestset "Eigh wrapper" begin
-            include("utility/eigh_wrapper.jl")
-        end
-        @time @safetestset "SVD wrapper" begin
-            include("utility/svd_wrapper.jl")
-        end
-        @time @safetestset "Symmetrization" begin
-            include("utility/symmetrization.jl")
-        end
-        @time @safetestset "Differentiable tmap" begin
-            include("utility/diff_maps.jl")
-        end
-        @time @safetestset "Norm-preserving retractions" begin
-            include("utility/retractions.jl")
-        end
-        @time @safetestset "Correlators" begin
-            include("utility/correlator.jl")
-        end
-    end
-    if GROUP == "ALL" || GROUP == "EXAMPLES"
-        @time @safetestset "Transverse field Ising model" begin
-            include("examples/tf_ising.jl")
-        end
-        @time @safetestset "Heisenberg model" begin
-            include("examples/heisenberg.jl")
-        end
-        @time @safetestset "J1-J2 model" begin
-            include("examples/j1j2_model.jl")
-        end
-        @time @safetestset "P-wave superconductor" begin
-            include("examples/pwave.jl")
-        end
-        @time @safetestset "U1-symmetric Bose-Hubbard model" begin
-            include("examples/bose_hubbard.jl")
-        end
-    end
+const init_code = quote
+    const fast_tests = $fast
 end
+
+ParallelTestRunner.runtests(PEPSKit, args; init_code)
diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl
index 1276ad2f8..cb054d56a 100644
--- a/test/timeevol/cluster_projectors.jl
+++ b/test/timeevol/cluster_projectors.jl
@@ -6,7 +6,123 @@ using Random
 import MPSKitModels: hubbard_space
 using PEPSKit: sdiag_pow, _cluster_truncate!, _flip_virtuals!, _next
 using MPSKit: GenericMPSTensor, MPSBondTensor
-include("cluster_tools.jl")
+
+# Utility setup
+# -------------
+function _contract_left(
+        M::GenericMPSTensor{S, 4}, sl::DiagonalTensorMap{T, S}
+    ) where {T <: Number, S <: ElementarySpace}
+    @assert !isdual(codomain(M, 1)) && !isdual(domain(M, 1))
+    M0 = twist(M, filter(ax -> isdual(space(M, ax)), 1:4))
+    @tensor sl1[e1; e0] := conj(M[w1; p n s e1]) * sl[w1; w0] * M0[w0; p n s e0]
+    return sl1
+end
+function _contract_left(
+        M::GenericMPSTensor{S, 4}, ::Nothing
+    ) where {S <: ElementarySpace}
+    @assert !isdual(domain(M, 1))
+    M0 = twist(M, filter(ax -> isdual(space(M, ax)), 1:4))
+    @tensor sl1[e1; e0] := conj(M[w; p n s e1]) * M0[w; p n s e0]
+    return sl1
+end
+
+function _contract_right(
+        M::GenericMPSTensor{S, 4}, sr::DiagonalTensorMap{T, S}
+    ) where {T <: Number, S <: ElementarySpace}
+    @assert !isdual(codomain(M, 1)) && !isdual(domain(M, 1))
+    M0 = twist(M, filter(ax -> !isdual(space(M, ax)), 2:5))
+    @tensor sr1[w0; w1] := M0[w0; p n s e0] * sr[e0; e1] * conj(M[w1; p n s e1])
+    return sr1
+end
+function _contract_right(
+        M::GenericMPSTensor{S, 4}, ::Nothing
+    ) where {S <: ElementarySpace}
+    @assert !isdual(codomain(M, 1))
+    M0 = twist(M, filter(ax -> !isdual(space(M, ax)), 2:5))
+    @tensor sr1[w0; w1] := M0[w0; p n s e] * conj(M[w1; p n s e])
+    return sr1
+end
+
+"""
+Verify the generalized left/right orthogonal condition
+"""
+function verify_cluster_orth(
+        Ms::Vector{T1}, wts::Vector{T2}
+    ) where {T1 <: GenericMPSTensor{<:ElementarySpace, 4}, T2 <: DiagonalTensorMap}
+    N = length(Ms)
+    @assert length(wts) == N - 1
+    lorths = fill(false, N - 1)
+    rorths = fill(false, N - 1)
+    # left orthogonal
+    for i in 1:(N - 1)
+        M, sl0 = Ms[i], wts[i]
+        sl1 = _contract_left(M, i == 1 ? nothing : wts[i - 1])
+        lorths[i] = (normalize(TensorMap(sl0)) ≈ normalize(sl1)) # sl0 is DiagonalTensorMap while sl1 is not
+    end
+    # right orthogonal
+    for i in 2:N
+        M, sr0 = Ms[i], wts[i - 1]
+        sr1 = _contract_right(M, i == N ? nothing : wts[i])
+        rorths[i - 1] = (normalize(TensorMap(sr0)) ≈ normalize(sr1))
+    end
+    return lorths, rorths
+end
+
+function inner_prod_cluster(
+        Ms1::Vector{T1}, Ms2::Vector{T2}
+    ) where {
+        T1 <: GenericMPSTensor{<:ElementarySpace, 4},
+        T2 <: GenericMPSTensor{<:ElementarySpace, 4},
+    }
+    N = length(Ms1)
+    @assert length(Ms2) == N
+    # physical spaces are assumed to be non-dual
+    @assert all(!isdual(space(t, 2)) for t in Ms1)
+    @assert all(!isdual(space(t, 2)) for t in Ms2)
+    # not the most efficient implementation
+    M1, M2 = Ms1[1], deepcopy(Ms2[1])
+    for ax in 1:4
+        isdual(space(M2, ax)) && twist!(M2, ax)
+    end
+    @tensor res[-1 -2] := conj(M1[1 2 3 4; -1]) * M2[1 2 3 4; -2]
+    for i in 2:(N - 1)
+        M1, M2 = Ms1[i], deepcopy(Ms2[i])
+        for ax in 2:4
+            isdual(space(M2, ax)) && twist!(M2, ax)
+        end
+        @tensor M[-1 -2; -3 -4] := conj(M1[-1 1 2 3; -3]) * M2[-2 1 2 3; -4]
+        @tensor res[-1 -2] := res[1 2] * M[1 2; -1 -2]
+    end
+    M1, M2 = Ms1[N], deepcopy(Ms2[N])
+    for ax in 2:5
+        isdual(space(M2, ax)) && twist!(M2, ax)
+    end
+    @tensor M[-1 -2] := conj(M1[-1 1 2 3; 4]) * M2[-2 1 2 3; 4]
+    return @tensor res[1 2] * M[1 2]
+end
+
+function fidelity_cluster(
+        Ms1::Vector{T1}, Ms2::Vector{T2}
+    ) where {
+        T1 <: GenericMPSTensor{<:ElementarySpace, 4},
+        T2 <: GenericMPSTensor{<:ElementarySpace, 4},
+    }
+    return abs2(inner_prod_cluster(Ms1, Ms2)) /
+        (inner_prod_cluster(Ms1, Ms1) * inner_prod_cluster(Ms2, Ms2))
+end
+
+function mpo_to_gate3(gs::Vector{T}) where {T <: AbstractTensorMap}
+    #= 
+    -4         -5        -6
+    ↓           ↓         ↓
+    g1 ←- 1 ←- g2 ←- 2 ←- g3
+    ↓           ↓         ↓
+    -1         -2        -3
+    =#
+    @assert length(gs) == 3
+    @tensor gate[-1 -2 -3; -4 -5 -6] := gs[1][-1 -4 1] * gs[2][1 -2 -5 2] * gs[3][2 -3 -6]
+    return gate
+end
 
 Vspaces = [
     (
diff --git a/test/timeevol/cluster_tools.jl b/test/timeevol/cluster_tools.jl
deleted file mode 100644
index efcdfc770..000000000
--- a/test/timeevol/cluster_tools.jl
+++ /dev/null
@@ -1,114 +0,0 @@
-function _contract_left(
-        M::GenericMPSTensor{S, 4}, sl::DiagonalTensorMap{T, S}
-    ) where {T <: Number, S <: ElementarySpace}
-    @assert !isdual(codomain(M, 1)) && !isdual(domain(M, 1))
-    M0 = twist(M, filter(ax -> isdual(space(M, ax)), 1:4))
-    @tensor sl1[e1; e0] := conj(M[w1; p n s e1]) * sl[w1; w0] * M0[w0; p n s e0]
-    return sl1
-end
-function _contract_left(
-        M::GenericMPSTensor{S, 4}, ::Nothing
-    ) where {S <: ElementarySpace}
-    @assert !isdual(domain(M, 1))
-    M0 = twist(M, filter(ax -> isdual(space(M, ax)), 1:4))
-    @tensor sl1[e1; e0] := conj(M[w; p n s e1]) * M0[w; p n s e0]
-    return sl1
-end
-
-function _contract_right(
-        M::GenericMPSTensor{S, 4}, sr::DiagonalTensorMap{T, S}
-    ) where {T <: Number, S <: ElementarySpace}
-    @assert !isdual(codomain(M, 1)) && !isdual(domain(M, 1))
-    M0 = twist(M, filter(ax -> !isdual(space(M, ax)), 2:5))
-    @tensor sr1[w0; w1] := M0[w0; p n s e0] * sr[e0; e1] * conj(M[w1; p n s e1])
-    return sr1
-end
-function _contract_right(
-        M::GenericMPSTensor{S, 4}, ::Nothing
-    ) where {S <: ElementarySpace}
-    @assert !isdual(codomain(M, 1))
-    M0 = twist(M, filter(ax -> !isdual(space(M, ax)), 2:5))
-    @tensor sr1[w0; w1] := M0[w0; p n s e] * conj(M[w1; p n s e])
-    return sr1
-end
-
-"""
-Verify the generalized left/right orthogonal condition
-"""
-function verify_cluster_orth(
-        Ms::Vector{T1}, wts::Vector{T2}
-    ) where {T1 <: GenericMPSTensor{<:ElementarySpace, 4}, T2 <: DiagonalTensorMap}
-    N = length(Ms)
-    @assert length(wts) == N - 1
-    lorths = fill(false, N - 1)
-    rorths = fill(false, N - 1)
-    # left orthogonal
-    for i in 1:(N - 1)
-        M, sl0 = Ms[i], wts[i]
-        sl1 = _contract_left(M, i == 1 ? nothing : wts[i - 1])
-        lorths[i] = (normalize(TensorMap(sl0)) ≈ normalize(sl1)) # sl0 is DiagonalTensorMap while sl1 is not
-    end
-    # right orthogonal
-    for i in 2:N
-        M, sr0 = Ms[i], wts[i - 1]
-        sr1 = _contract_right(M, i == N ? nothing : wts[i])
-        rorths[i - 1] = (normalize(TensorMap(sr0)) ≈ normalize(sr1))
-    end
-    return lorths, rorths
-end
-
-function inner_prod_cluster(
-        Ms1::Vector{T1}, Ms2::Vector{T2}
-    ) where {
-        T1 <: GenericMPSTensor{<:ElementarySpace, 4},
-        T2 <: GenericMPSTensor{<:ElementarySpace, 4},
-    }
-    N = length(Ms1)
-    @assert length(Ms2) == N
-    # physical spaces are assumed to be non-dual
-    @assert all(!isdual(space(t, 2)) for t in Ms1)
-    @assert all(!isdual(space(t, 2)) for t in Ms2)
-    # not the most efficient implementation
-    M1, M2 = Ms1[1], deepcopy(Ms2[1])
-    for ax in 1:4
-        isdual(space(M2, ax)) && twist!(M2, ax)
-    end
-    @tensor res[-1 -2] := conj(M1[1 2 3 4; -1]) * M2[1 2 3 4; -2]
-    for i in 2:(N - 1)
-        M1, M2 = Ms1[i], deepcopy(Ms2[i])
-        for ax in 2:4
-            isdual(space(M2, ax)) && twist!(M2, ax)
-        end
-        @tensor M[-1 -2; -3 -4] := conj(M1[-1 1 2 3; -3]) * M2[-2 1 2 3; -4]
-        @tensor res[-1 -2] := res[1 2] * M[1 2; -1 -2]
-    end
-    M1, M2 = Ms1[N], deepcopy(Ms2[N])
-    for ax in 2:5
-        isdual(space(M2, ax)) && twist!(M2, ax)
-    end
-    @tensor M[-1 -2] := conj(M1[-1 1 2 3; 4]) * M2[-2 1 2 3; 4]
-    return @tensor res[1 2] * M[1 2]
-end
-
-function fidelity_cluster(
-        Ms1::Vector{T1}, Ms2::Vector{T2}
-    ) where {
-        T1 <: GenericMPSTensor{<:ElementarySpace, 4},
-        T2 <: GenericMPSTensor{<:ElementarySpace, 4},
-    }
-    return abs2(inner_prod_cluster(Ms1, Ms2)) /
-        (inner_prod_cluster(Ms1, Ms1) * inner_prod_cluster(Ms2, Ms2))
-end
-
-function mpo_to_gate3(gs::Vector{T}) where {T <: AbstractTensorMap}
-    #= 
-    -4         -5        -6
-    ↓           ↓         ↓
-    g1 ←- 1 ←- g2 ←- 2 ←- g3
-    ↓           ↓         ↓
-    -1         -2        -3
-    =#
-    @assert length(gs) == 3
-    @tensor gate[-1 -2 -3; -4 -5 -6] := gs[1][-1 -4 1] * gs[2][1 -2 -5 2] * gs[3][2 -3 -6]
-    return gate
-end