diff --git a/Project.toml b/Project.toml
index 4e9564f1..74c7a8b4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -15,6 +15,7 @@ OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
@@ -28,6 +29,7 @@ MatrixAlgebraKit = "0.6.1"
 OptimKit = "0.4.0"
 QuadGK = "2.11.2"
 SpecialFunctions = "2.6.1"
+StableRNGs = "1.0.4"
 TensorKit = "0.16.2"
 Zygote = "0.7.7"
 julia = "1.11"
diff --git a/README.md b/README.md
index ca0587fb..20295c22 100644
--- a/README.md
+++ b/README.md
@@ -48,6 +48,7 @@ The following schemes are currently implemented:
 
 **2D honeycomb CTM methods**
 - c3vCTM_honeycomb (c3v symmetric CTM on the honeycomb lattice)
+- CTM_honeycomb (CTM on the honeycomb lattice)
 
 **2D Impurity Methods**
 - ImpurityTRG (Expectation value calculation via TRG)
diff --git a/docs/src/assets/tnrkit.bib b/docs/src/assets/tnrkit.bib
index 422916e0..97a09caf 100644
--- a/docs/src/assets/tnrkit.bib
+++ b/docs/src/assets/tnrkit.bib
@@ -224,6 +224,21 @@ @article{lukin2023
   url       = {https://link.aps.org/doi/10.1103/PhysRevB.107.054424}
 }
 
+@article{nyckees2023,
+  title = {Critical line of the triangular Ising antiferromagnet in a field from a ${C}_{3}$-symmetric corner transfer matrix algorithm},
+  author = {Nyckees, Samuel and Rufino, Afonso and Mila, Fr\'ed\'eric and Colbois, Jeanne},
+  journal = {Phys. Rev. E},
+  volume = {108},
+  issue = {6},
+  pages = {064132},
+  numpages = {14},
+  year = {2023},
+  month = {Dec},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevE.108.064132},
+  url = {https://link.aps.org/doi/10.1103/PhysRevE.108.064132}
+}
+
 @article{homma2024a,
   title   = {Nuclear Norm Regularized Loop Optimization for Tensor Network},
   author  = {Homma, Kenji},
diff --git a/docs/src/index.md b/docs/src/index.md
index a4bec77f..0d77b1f2 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -33,6 +33,7 @@ Many common TNR schemes have already been implemented:
 
 **2D honeycomb CTM methods**
 * [`c3vCTM_honeycomb`](@ref) (c3v symmetric CTM on the honeycomb lattice)
+* [`CTM_honeycomb`](@ref) (CTM on the honeycomb lattice)
 
 **Impurity Methods**
 * [`ImpurityTRG`](@ref)
diff --git a/src/TNRKit.jl b/src/TNRKit.jl
index 0c32cbc0..cc24a258 100644
--- a/src/TNRKit.jl
+++ b/src/TNRKit.jl
@@ -36,6 +36,8 @@ include("schemes/ctm/sublattice_ctm.jl")
 include("schemes/ctm/triangular.jl")
 include("schemes/ctm/ctm_triangular.jl")
 include("schemes/ctm/c6vctm_triangular.jl")
+include("schemes/ctm/honeycomb.jl")
+include("schemes/ctm/ctm_honeycomb.jl")
 include("schemes/ctm/c3vctm_honeycomb.jl")
 
 # Impurity methods
@@ -69,6 +71,7 @@ export lnz
 export c6vCTM_triangular
 export CTM_triangular
 export c3vCTM_honeycomb
+export CTM_honeycomb
 
 export ImpurityTRG
 export ImpurityHOTRG
diff --git a/src/schemes/ctm/c3vctm_honeycomb.jl b/src/schemes/ctm/c3vctm_honeycomb.jl
index 7806ac3f..6793fc03 100644
--- a/src/schemes/ctm/c3vctm_honeycomb.jl
+++ b/src/schemes/ctm/c3vctm_honeycomb.jl
@@ -1,94 +1,3 @@
-"""
-$(TYPEDEF)
-
-Corner Transfer Matrix Renormalization Group for the honeycomb lattice
-
-### Constructors
-    $(FUNCTIONNAME)(T)
-    $(FUNCTIONNAME)(T, [, symmetrize=false])
-
-```
-     (120°)
-        ╲ 
-         ╲ 
-          ╲ 
-           T -----(0°)
-           ╱
-          ╱
-         ╱
-      (240°)
-```
-
-CTM can be called with a (2, 1) tensor, where the directions are (240°, 0°, 120°) clockwise with respect to the positive x-axis.
-In the flipped arrow convention, the arrows point from (120°) to (240°, 0°).
-or with a (0,3) tensor (120°, 0°, 240°) where all arrows point inward (unflipped arrow convention).
-The keyword argument symmetrize makes the tensor C6v symmetric when set to true. If symmetrize = false, it checks the symmetry explicitly.
-
-### Running the algorithm
-    run!(::CTM, trunc::TruncationStrategy, stop::Stopcrit[, finalize_beginning=true, verbosity=1])
-
-!!! info "verbosity levels"
-    - 0: No output
-    - 1: Print information at start and end of the algorithm
-    - 2: Print information at each step
-
-### Fields
-
-$(TYPEDFIELDS)
-
-### References
-* [Lukin et al. Phys. Rev. B 107.054424 (2023)](@cite lukin2023)
-"""
-mutable struct c3vCTM_honeycomb{A, S}
-    T::TensorMap{A, S, 0, 3}
-    C::TensorMap{A, S, 1, 1}
-    L::TensorMap{A, S, 2, 1}
-    R::TensorMap{A, S, 2, 1}
-
-    function c3vCTM_honeycomb(T::TensorMap{A, S, 0, 3}) where {A, S}
-        C, Ea, Eb = c3vCTM_honeycomb_init(T)
-
-        if BraidingStyle(sectortype(T)) != Bosonic()
-            @warn "$(summary(BraidingStyle(sectortype(T)))) braiding style is not supported for c6vCTM"
-        end
-        return new{A, S}(T, C, Ea, Eb)
-    end
-end
-
-function c3vCTM_honeycomb(T_flipped::TensorMap{A, S, 2, 1}; symmetrize = false) where {A, S}
-    T_unflipped = permute(flip(T_flipped, [1 2]; inv = true), ((), (3, 2, 1)))
-    if symmetrize
-        T_unflipped = symmetrize_C6v_honeycomb(T_unflipped)
-    else
-        @assert norm(T_unflipped - rotl120_pf_honeycomb(T_unflipped)) < 1.0e-14 "Tensor is not C6 symmetric. Error = $(norm(T_unflipped - rotl60_pf(T_unflipped)))"
-    end
-
-    return c3vCTM_honeycomb(T_unflipped)
-end
-
-function c3vCTM_honeycomb_init(T::TensorMap{A, S, 0, 3}) where {A, S}
-    S_type = scalartype(T)
-    Vp = space(T)[1]'
-    C = ones(S_type, oneunit(Vp) ← oneunit(Vp))
-    L = ones(S_type, oneunit(Vp) ⊗ Vp ← oneunit(Vp))
-    R = ones(S_type, oneunit(Vp) ⊗ Vp ← oneunit(Vp))
-    return C, L, R
-end
-
-# Functions to permute (flipped and unflipped) tensors under 60 degree rotation
-function rotl120_pf_honeycomb(T::TensorMap{A, S, 2, 1}) where {A, S}
-    return permute(T, ((3, 1), (2,)))
-    return permute(T, ((4, 1, 2), (5, 6, 3)))
-end
-
-function rotl120_pf_honeycomb(T::TensorMap{A, S, 0, 3}) where {A, S}
-    return permute(T, ((), (2, 3, 1)))
-end
-
-function symmetrize_C6v_honeycomb(T_unflipped)
-    return (T_unflipped + rotl120_pf_honeycomb(T_unflipped) + rotl120_pf_honeycomb(rotl120_pf_honeycomb(T_unflipped))) / 3
-end
-
 # Based on
 # https://arxiv.org/abs/2209.03428
 
diff --git a/src/schemes/ctm/ctm_honeycomb.jl b/src/schemes/ctm/ctm_honeycomb.jl
new file mode 100644
index 00000000..5d6f56c3
--- /dev/null
+++ b/src/schemes/ctm/ctm_honeycomb.jl
@@ -0,0 +1,134 @@
+function run!(
+        scheme::CTM_honeycomb, trunc::TruncationStrategy, criterion::stopcrit;
+        verbosity = 1
+    )
+    LoggingExtras.withlevel(; verbosity) do
+        @infov 1 "Starting simulation\n"
+        steps = 0
+        crit = true
+        ε = Inf
+        Ss_prev = [DiagonalTensorMap(id(domain(scheme.C[dir]))) for dir in 1:3]
+        t = @elapsed while crit
+            @infov 2 "Step $(steps + 1), ε = $(ε)"
+
+            Ss = step!(scheme, trunc)
+            ε = error_measure(Ss, Ss_prev)
+            Ss_prev = Ss
+            steps += 1
+            crit = criterion(steps, ε)
+        end
+
+        @infov 1 "Simulation finished\n $(stopping_info(criterion, steps, ε))\n Elapsed time: $(t)s\n Iterations: $steps"
+    end
+    return lnz(scheme)
+end
+
+function step!(scheme::CTM_honeycomb, trunc)
+    # take two half-steps to compensate for the swapping of Ta and Tb in one half-step.
+    half_step!(scheme, trunc)
+    Ss = half_step!(scheme, trunc)
+    return Ss
+end
+
+function half_step!(scheme::CTM_honeycomb, trunc)
+    projectors, Ss = calculate_projectors(scheme, trunc)
+    renormalize_edges!(scheme, projectors)
+    return Ss
+end
+
+function build_corner_matrices(scheme::CTM_honeycomb)
+    @tensor opt = true σ₁[χNW1 DaSW; χNE2 DbSE] := flip(scheme.A, [2 3])[DaNW DaE DaSW] * scheme.B[DbSE DaE DbNE] *
+        scheme.Tb[3][χNW1 DaNW; χNW2] * scheme.Ta[1][χNE1 DbNE; χNE2] *
+        scheme.C[1][χNW2; χNE1]
+    @tensor opt = true σ₂[χNE1 DaNW; χS2 DbW] := flip(scheme.A, [1 3])[DaNW DaE DaSW] * scheme.B[DbSE DbW DaSW] *
+        scheme.Tb[1][χNE1 DaE; χNE2] * scheme.Ta[2][χS1 DbSE; χS2] *
+        scheme.C[2][χNE2; χS1]
+    @tensor opt = true σ₃[χS1 DaE; χNW2 DbNE] := flip(scheme.A, [1 2])[DaNW DaE DaSW] * scheme.B[DaNW DbW DbNE] *
+        scheme.Tb[2][χS1 DaSW; χS2] * scheme.Ta[3][χNW1 DbW; χNW2] *
+        scheme.C[3][χS2; χNW1]
+    return [σ₁, σ₂, σ₃]
+end
+
+function calculate_projector_simultaneous(σs, trunc, dir)
+    σsshift = circshift(σs, 1 - dir)
+    mat = prod(σsshift)
+    U, S, Vᴴ = svd_trunc(mat; trunc)
+    PL = prod(σsshift[2:end]) * Vᴴ' * pseudopow(S, -1 / 2)
+    PR = pseudopow(S, -1 / 2) * U' * σsshift[1]
+    return PL, PR, S
+end
+
+function calculate_projectors(scheme::CTM_honeycomb, trunc)
+    σs = build_corner_matrices(scheme)
+    projectors_simul = [calculate_projector_simultaneous(σs, trunc, dir) for dir in 1:3]
+    projectors = [projectors_simul[dir][1:2] for dir in 1:3]
+    Ss = [projectors_simul[dir][3] for dir in 1:3]
+    for dir in 1:3
+        C′ = projectors[mod1(dir - 1, 3)][2] * σs[dir] * projectors[dir][1]
+        scheme.C[dir] = C′ / norm(C′)
+    end
+    return projectors, Ss
+end
+
+function renormalize_edges!(scheme::CTM_honeycomb, projectors)
+    for dir in 1:3
+        @tensor opt = true Ta′[χout DSW; χin] := projectors[dir][2][χout; χNE1 DNW] * scheme.Tb[dir][χNE1 DE; χin] * flip(rotl120_pf_honeycomb(scheme.A, dir - 1), [1 3])[DNW DE DSW]
+        @tensor opt = true Tb′[χout DW; χin] := scheme.Ta[dir][χout DNE; χNE] * flip(rotl120_pf_honeycomb(scheme.B, dir - 1), 2)[DSE DW DNE] * projectors[dir][1][χNE DSE; χin]
+        scheme.Ta[dir] = Ta′ / norm(Ta′)
+        scheme.Tb[dir] = Tb′ / norm(Tb′)
+    end
+    return
+end
+
+function error_measure(Ss::Vector{T}, Ss_prev::Vector{T}) where {E, U, T <: AbstractTensorMap{E, U}}
+    ϵs = [(space(S) == space(S_prev)) ? norm(S - S_prev) : Inf for (S, S_prev) in zip(Ss, Ss_prev)]
+    return maximum(ϵs)
+end
+
+function lnz(scheme::CTM_honeycomb)
+    return real(log(network_value(scheme)))
+end
+
+function network_value(scheme::CTM_honeycomb)
+    nw_full_large = _contract_site_large(scheme)
+    nw_full = _contract_site(scheme)
+    nw_corners = _contract_corners(scheme)
+    return (nw_full_large * nw_corners / nw_full^2)^(1 / 6)
+end
+
+function _contract_site_large(scheme::CTM_honeycomb)
+    return @tensor opt = true scheme.C[1][χNW5; χNE1] *
+        scheme.Ta[1][χNE1 D2; χNE2] * scheme.Tb[1][χNE2 D7; χNE3] *
+        scheme.Ta[1][χNE3 D14; χNE4] * scheme.Tb[1][χNE4 D23; χNE5] *
+        scheme.C[2][χNE5; χS1] *
+        scheme.Ta[2][χS1 D33; χS2] * scheme.Tb[2][χS2 D32; χS3] *
+        scheme.Ta[2][χS3 D31; χS4] * scheme.Tb[2][χS4 D30; χS5] *
+        scheme.C[3][χS5; χNW1] *
+        scheme.Ta[3][χNW1 D22; χNW2] * scheme.Tb[3][χNW2 D13; χNW3] *
+        scheme.Ta[3][χNW3 D6; χNW4] * scheme.Tb[3][χNW4 D1; χNW5] *
+        scheme.A[D1 D3 D4] * flip(scheme.B, [1 2])[D5 D3 D2] *
+        scheme.A[D5 D7 D9] * flip(scheme.B, [1 2 3])[D12 D10 D9] *
+        scheme.A[D8 D10 D11] * flip(scheme.B, [1 3])[D8 D6 D4] *
+        scheme.A[D13 D15 D17] * flip(scheme.B, [1 2 3])[D18 D15 D11] *
+        scheme.A[D18 D21 D25] * flip(scheme.B, [2 3])[D31 D28 D25] *
+        scheme.A[D24 D28 D30] * flip(scheme.B, [1 3])[D24 D22 D17] *
+        scheme.A[D12 D16 D19] * flip(scheme.B, [1 2])[D20 D16 D14] *
+        scheme.A[D20 D23 D27] * flip(scheme.B, [2 3])[D33 D29 D27] *
+        scheme.A[D26 D29 D32] * flip(scheme.B, [1 2 3])[D26 D21 D19]
+end
+
+function _contract_site(scheme::CTM_honeycomb)
+    return @tensor opt = true scheme.C[1][χNW2; χNE1] *
+        scheme.Ta[1][χNE1 D1; χNEC] * scheme.Tb[1][χNEC D2; χNE2] *
+        scheme.C[2][χNE2; χS1] *
+        scheme.Ta[2][χS1 D3; χSC] * scheme.Tb[2][χSC D4; χS2] *
+        scheme.C[3][χS2; χNW1] *
+        scheme.Ta[3][χNW1 D5; χNWC] * scheme.Tb[3][χNWC D6; χNW2] *
+        scheme.A[D6 D7 D12] * flip(scheme.B, [1 2])[D8 D7 D1] *
+        scheme.A[D8 D2 D9] * flip(scheme.B, [2 3])[D3 D10 D9] *
+        scheme.A[D11 D10 D4] * flip(scheme.B, [1 3])[D11 D5 D12]
+end
+
+function _contract_corners(scheme::CTM_honeycomb)
+    return @tensor scheme.C[1][1; 2] * scheme.C[2][2; 3] * scheme.C[3][3; 1]
+end
diff --git a/src/schemes/ctm/honeycomb.jl b/src/schemes/ctm/honeycomb.jl
new file mode 100644
index 00000000..ed41a74a
--- /dev/null
+++ b/src/schemes/ctm/honeycomb.jl
@@ -0,0 +1,168 @@
+"""
+$(TYPEDEF)
+
+Corner Transfer Matrix Renormalization Group for the honeycomb lattice
+
+### Constructors
+    $(FUNCTIONNAME)(T)
+    $(FUNCTIONNAME)(T, [, symmetrize=false])
+
+```
+     (120°)
+        ╲ 
+         ╲ 
+          ╲ 
+           T -----(0°)
+           ╱
+          ╱
+         ╱
+      (240°)
+```
+
+CTM can be called with a (2, 1) tensor, where the directions are (240°, 0°, 120°) clockwise with respect to the positive x-axis.
+In the flipped arrow convention, the arrows point from (120°) to (240°, 0°).
+or with a (0,3) tensor (120°, 0°, 240°) where all arrows point inward (unflipped arrow convention).
+The keyword argument symmetrize makes the tensor C6v symmetric when set to true. If symmetrize = false, it checks the symmetry explicitly.
+
+### Running the algorithm
+    run!(::CTM, trunc::TruncationStrategy, stop::Stopcrit[, finalize_beginning=true, verbosity=1])
+
+!!! info "verbosity levels"
+    - 0: No output
+    - 1: Print information at start and end of the algorithm
+    - 2: Print information at each step
+
+### Fields
+
+$(TYPEDFIELDS)
+
+### References
+* [Lukin et al. Phys. Rev. B 107.054424 (2023)](@cite lukin2023)
+"""
+mutable struct c3vCTM_honeycomb{A, S}
+    T::TensorMap{A, S, 0, 3}
+    C::TensorMap{A, S, 1, 1}
+    L::TensorMap{A, S, 2, 1}
+    R::TensorMap{A, S, 2, 1}
+
+    function c3vCTM_honeycomb(T::TensorMap{A, S, 0, 3}) where {A, S}
+        C, Ea, Eb = c3vCTM_honeycomb_init(T)
+
+        if BraidingStyle(sectortype(T)) != Bosonic()
+            @warn "$(summary(BraidingStyle(sectortype(T)))) braiding style is not supported for c6vCTM"
+        end
+        return new{A, S}(T, C, Ea, Eb)
+    end
+end
+
+function c3vCTM_honeycomb(T_flipped::TensorMap{A, S, 2, 1}; symmetrize = false) where {A, S}
+    T_unflipped = permute(flip(T_flipped, [1 2]; inv = true), ((), (3, 2, 1)))
+    if symmetrize
+        T_unflipped = symmetrize_C3_honeycomb(T_unflipped)
+    else
+        !is_C3_symmetric(T_unflipped) && throw(ArgumentError("Tensor is not C3 symmetric"))
+    end
+    return c3vCTM_honeycomb(T_unflipped)
+end
+
+function c3vCTM_honeycomb_init(T::TensorMap{A, S, 0, 3}) where {A, S}
+    S_type = scalartype(T)
+    Vp = space(T)[1]'
+    C = ones(S_type, oneunit(Vp) ← oneunit(Vp))
+    L = ones(S_type, oneunit(Vp) ⊗ Vp ← oneunit(Vp))
+    R = ones(S_type, oneunit(Vp) ⊗ Vp ← oneunit(Vp))
+    return C, L, R
+end
+
+# Functions to permute (unflipped) tensors under 120 degree rotation
+function rotl120_pf_honeycomb(T::TensorMap{A, S, 0, 3}) where {A, S}
+    return permute(T, ((), (2, 3, 1)))
+end
+
+function rotl120_pf_honeycomb(T::TensorMap{A, S, 0, 3}, i::Int) where {A, S}
+    @assert i >= 0 "Negative rotation not supported"
+    if i == 0
+        return T
+    end
+    return rotl120_pf_honeycomb(rotl120_pf_honeycomb(T), i - 1)
+end
+
+function symmetrize_C3_honeycomb(T_unflipped::TensorMap{E, S, 0, 3}) where {E, S}
+    return (T_unflipped + rotl120_pf_honeycomb(T_unflipped) + rotl120_pf_honeycomb(rotl120_pf_honeycomb(T_unflipped))) / 3
+end
+
+function is_C3_symmetric(T_unflipped::TensorMap{E, S, 0, 3}) where {E, S}
+    return space(T_unflipped) == space(rotl120_pf_honeycomb(T_unflipped)) && norm(T_unflipped - rotl120_pf_honeycomb(T_unflipped)) < 1.0e-14
+end
+
+
+"""
+$(TYPEDEF)
+
+Corner Transfer Matrix Renormalization Group for the honeycomb lattice
+
+### Constructors
+    $(FUNCTIONNAME)(T)
+    $(FUNCTIONNAME)(T, [, symmetrize=false])
+
+```
+     (120°)
+        ╲ 
+         ╲ 
+          ╲ 
+           T -----(0°)
+           ╱
+          ╱
+         ╱
+      (240°)
+```
+
+CTM can be called with a (2, 1) tensor, where the directions are (240°, 0°, 120°) clockwise with respect to the positive x-axis.
+In the flipped arrow convention, the arrows point from (120°) to (240°, 0°).
+or with a (0,3) tensor (120°, 0°, 240°) where all arrows point inward (unflipped arrow convention).
+
+### Running the algorithm
+    run!(::CTM, trunc::TruncationStrategy, stop::Stopcrit[, finalize_beginning=true, verbosity=1])
+
+!!! info "verbosity levels"
+    - 0: No output
+    - 1: Print information at start and end of the algorithm
+    - 2: Print information at each step
+
+### Fields
+
+$(TYPEDFIELDS)
+
+### References
+* [Nyckees et al. Phys. Rev. E 108.064132 (2023)](@cite nyckees2023)
+"""
+mutable struct CTM_honeycomb{E, S}
+    A::TensorMap{E, S, 0, 3}
+    B::TensorMap{E, S, 0, 3}
+    C::Vector{TensorMap{E, S, 1, 1}}
+    Ta::Vector{TensorMap{E, S, 2, 1}}
+    Tb::Vector{TensorMap{E, S, 2, 1}}
+
+    function CTM_honeycomb(A::TensorMap{E, S, 0, 3}; B::TensorMap{E, S, 0, 3} = A) where {E, S}
+        C, Ta, Tb = CTM_honeycomb_init(A; B)
+        if BraidingStyle(sectortype(A)) != Bosonic()
+            @warn "$(summary(BraidingStyle(sectortype(A)))) braiding style is not supported for c6vCTM"
+        end
+        return new{E, S}(A, B, C, Ta, Tb)
+    end
+end
+
+function CTM_honeycomb(A_flipped::TensorMap{E, S, 2, 1}; B::TensorMap{E, S, 2, 1} = A_flipped) where {E, S}
+    A_unflipped = permute(flip(A_flipped, [1 2]; inv = true), ((), (3, 2, 1)))
+    B_unflipped = permute(flip(B, [1 2]; inv = true), ((), (3, 2, 1)))
+    return CTM_honeycomb(A_unflipped; B = B_unflipped)
+end
+
+function CTM_honeycomb_init(A::TensorMap{E, S, 0, 3}; B::TensorMap{E, S, 0, 3} = A) where {E, S}
+    S_type = scalartype(A)
+    Vp = space(A)[1]'
+    C = fill(ones(S_type, oneunit(Vp) ← oneunit(Vp)), 3)
+    Ta = [ones(S_type, oneunit(Vp) ⊗ space(B)[mod1(dir - 1, 3)]' ← oneunit(Vp)) for dir in 1:3]
+    Tb = [ones(S_type, oneunit(Vp) ⊗ space(A)[mod1(dir + 1, 3)]' ← oneunit(Vp)) for dir in 1:3]
+    return C, Ta, Tb
+end
diff --git a/test/schemes_honeycomb.jl b/test/schemes_honeycomb.jl
index e8643a70..54135659 100644
--- a/test/schemes_honeycomb.jl
+++ b/test/schemes_honeycomb.jl
@@ -1,12 +1,72 @@
-# c6CTM
-@testset "c6CTM - Ising Model" begin
+using StableRNGs
+@testset "Honeycomb schemes - Ising Model" begin
     for sym in [Trivial, Z2Irrep]
-        T_flipped = classical_ising_honeycomb(sym, ising_βc_honeycomb; T = ComplexF64)
-        println(T_flipped.space)
-        scheme = c3vCTM_honeycomb(T_flipped)
-        lz = run!(scheme, truncrank(20), convcrit(1.0e-4, (steps, data) -> data) & maxiter(300); verbosity = 1)
+        for alg in [:CTM_honeycomb, :c3vCTM_honeycomb]
+            T_flipped = classical_ising_honeycomb(sym, ising_βc_honeycomb; T = ComplexF64)
+            scheme = eval(alg)(T_flipped)
+            lz = run!(scheme, truncrank(20), convcrit(1.0e-4, (steps, data) -> data) & maxiter(300); verbosity = 1)
 
-        fs = lz * -1 / ising_βc_honeycomb
-        @test fs ≈ f_onsager_honeycomb rtol = 1.0e-2
+            fs = lz * -1 / ising_βc_honeycomb
+            @test fs ≈ f_onsager_honeycomb rtol = 1.0e-2
+        end
     end
 end
+
+# Test honeycomb CTM by converting it to CTM on a square lattice
+@testset "Honeycomb CTM C3 - Random Model" begin
+    rng = StableRNG(1234)
+    for sym in [Trivial, Z2Irrep]
+        for alg in [:CTM_honeycomb, :c3vCTM_honeycomb]
+            A = zeros(ComplexF64, ℂ^2 ⊗ ℂ^2, ℂ^2)
+            A.data .= rand(rng, length(A.data))
+            A /= norm(A)
+
+            if alg == :c3vCTM_honeycomb
+                scheme = eval(alg)(A; symmetrize = true)
+            else
+                scheme = eval(alg)(A)
+            end
+            lz_honeycomb = run!(scheme, truncrank(20), convcrit(1.0e-4, (steps, data) -> data) & maxiter(300); verbosity = 1)
+
+            @tensor pf_square[-4 -3; -1 -2] := A[-1 -2 1] * flip(A, [1 2 3])[-3 -4 1]
+            scheme_square = CTM(pf_square)
+            lz_square = run!(scheme, truncrank(20), convcrit(1.0e-4, (steps, data) -> data) & maxiter(300); verbosity = 1)
+
+            @test lz_square ≈ lz_honeycomb rtol = 1.0e-3
+        end
+    end
+end
+
+# Test CTM_honeycomb for A ≠ B by converting it to CTM on a square lattice
+@testset "Honeycomb CTM - Random Model" begin
+    rng = StableRNG(1234)
+    for sym in [Trivial, Z2Irrep]
+        A = zeros(ComplexF64, ℂ^2 ⊗ ℂ^3, ℂ^4)
+        B = zeros(ComplexF64, ℂ^2 ⊗ ℂ^3, ℂ^4)
+        A.data .= rand(rng, length(A.data))
+        B.data .= rand(rng, length(B.data))
+        A /= norm(A)
+        B /= norm(B)
+
+        @test_throws ArgumentError c3vCTM_honeycomb(A)
+        scheme = CTM_honeycomb(A; B)
+        lz_honeycomb = run!(scheme, truncrank(40), convcrit(-Inf, (steps, data) -> data) & maxiter(500); verbosity = 1)
+
+        @tensor pf_square[-4 -3; -1 -2] := A[-1 -2 1] * flip(B, [1 2 3])[-3 -4 1]
+        scheme_square = CTM(pf_square)
+        lz_square = run!(scheme, truncrank(40), convcrit(1.0e-14, (steps, data) -> data) & maxiter(500); verbosity = 1)
+
+        @test lz_square ≈ lz_honeycomb rtol = 1.0e-2
+    end
+end
+
+@testset "Rotations for honeycomb lattice" begin
+    rng = StableRNG(1234)
+    A_flipped = zeros(ComplexF64, ℂ^2 ⊗ ℂ^2, ℂ^2)
+    A = permute(flip(A_flipped, [1 2]; inv = true), ((), (3, 2, 1)))
+    A.data .= rand(rng, length(A.data))
+
+    @test TNRKit.rotl120_pf_honeycomb(A, 3) ≈ A
+    @test TNRKit.rotl120_pf_honeycomb(A, 1) ≈ TNRKit.rotl120_pf_honeycomb(A)
+    @test TNRKit.is_C3_symmetric(TNRKit.symmetrize_C3_honeycomb(A))
+end