From cb6b2a3290c9f43d5ab65d80214aaae710646286 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Fri, 26 Sep 2025 08:50:01 +0200
Subject: [PATCH 01/24] Add `initialize_environment`

---
 src/PEPSKit.jl                         |  3 ++
 src/algorithms/ctmrg/initialization.jl | 52 ++++++++++++++++++++++++++
 src/networks/infinitesquarenetwork.jl  |  1 +
 3 files changed, 56 insertions(+)
 create mode 100644 src/algorithms/ctmrg/initialization.jl

diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 9e1c24b74..1ae7f42a7 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -64,6 +64,7 @@ include("algorithms/ctmrg/projectors.jl")
 include("algorithms/ctmrg/simultaneous.jl")
 include("algorithms/ctmrg/sequential.jl")
 include("algorithms/ctmrg/gaugefix.jl")
+include("algorithms/ctmrg/initialization.jl")
 
 include("algorithms/truncation/truncationschemes.jl")
 include("algorithms/truncation/fullenv_truncation.jl")
@@ -87,6 +88,8 @@ using .Defaults: set_scheduler!
 export set_scheduler!
 export SVDAdjoint, FullSVDReverseRule, IterSVD
 export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG
+export initialize_environment,
+    RandomInitialization, ProductStateInitialization, ApplicationInitialization
 export FixedSpaceTruncation, SiteDependentTruncation
 export HalfInfiniteProjector, FullInfiniteProjector
 export LocalOperator, physicalspace
diff --git a/src/algorithms/ctmrg/initialization.jl b/src/algorithms/ctmrg/initialization.jl
new file mode 100644
index 000000000..cb5a5eebf
--- /dev/null
+++ b/src/algorithms/ctmrg/initialization.jl
@@ -0,0 +1,52 @@
+abstract type InitializationStyle end
+struct ProductStateInitialization <: InitializationStyle end
+struct RandomInitialization <: InitializationStyle end
+struct ApplicationInitialization <: InitializationStyle end
+
+function initialize_environment(
+        elt::Type{<:Number},
+        n::InfiniteSquareNetwork,
+        ::RandomInitialization,
+        init_spec::ElementarySpace = oneunit(spacetype(n)), # TODO: non-uniform space specification?
+    )
+    return CTMRGEnv(randn, elt, n, init_spec)
+end
+
+function initialize_environment(
+        elt::Type{<:Number},
+        n::InfiniteSquareNetwork,
+        ::ProductStateInitialization,
+        init_spec::ElementarySpace = oneunit(spacetype(n)), # TODO: non-uniform space specification?
+    )
+    i = one(sectortype(init_spec))
+    env = CTMRGEnv(ones, elt, n, init_spec)
+    for (dir, r, c) in Iterators.product(axes(env)...)
+        @assert i in blocksectors(env.corners[dir, r, c])
+        block(env.corners[dir, r, c], i)[1, 1] = 1
+    end
+    return env
+end
+
+function initialize_environment(
+        elt::Type{<:Number},
+        n::InfiniteSquareNetwork,
+        ::ApplicationInitialization,
+        init_spec::TruncationScheme;
+        boundary_alg = (;
+            alg = :sequential, tol = 1.0e-5, maxiter = 10, verbosity = -1, trscheme = init_spec,
+        )
+    )
+    env = initialize_environment(elt, n, ProductStateInitialization())
+    env, = leading_boundary(env, n; boundary_alg...)
+    return env
+end
+
+function initialize_environment(n::InfiniteSquareNetwork, args...; kwargs...)
+    return initialize_environment(ComplexF64, n, args...; kwargs...)
+end
+function initialize_environment(A::Union{InfinitePEPS, InfinitePartitionFunction}, args...; kwargs...)
+    return initialize_environment(ComplexF64, A, args...; kwargs...)
+end
+function initialize_environment(elt::Type{<:Number}, A::Union{InfinitePEPS, InfinitePartitionFunction}, args...; kwargs...)
+    return initialize_environment(elt, InfiniteSquareNetwork(A), args...; kwargs...)
+end
diff --git a/src/networks/infinitesquarenetwork.jl b/src/networks/infinitesquarenetwork.jl
index f003bd095..acb03c036 100644
--- a/src/networks/infinitesquarenetwork.jl
+++ b/src/networks/infinitesquarenetwork.jl
@@ -55,6 +55,7 @@ end
 
 ## Spaces
 
+TensorKit.spacetype(::Type{T}) where {T <: InfiniteSquareNetwork} = spacetype(eltype(T))
 virtualspace(n::InfiniteSquareNetwork, r::Int, c::Int, dir) = virtualspace(n[r, c], dir)
 
 ## Vector interface

From c1216f49098f3f91cc847ee0cbf3ec02a60cafc5 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Fri, 26 Sep 2025 11:19:39 +0200
Subject: [PATCH 02/24] Add test

---
 src/PEPSKit.jl                |  2 +-
 src/algorithms/ctmrg/ctmrg.jl |  2 ++
 test/ctmrg/initialization.jl  | 51 +++++++++++++++++++++++++++++++++++
 test/runtests.jl              |  3 +++
 4 files changed, 57 insertions(+), 1 deletion(-)
 create mode 100644 test/ctmrg/initialization.jl

diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 1ae7f42a7..a1363d541 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -2,7 +2,7 @@ module PEPSKit
 
 using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf
 using Compat
-using Accessors: @set, @reset
+using Accessors: @set, @reset, @insert
 using VectorInterface
 import VectorInterface as VI
 
diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl
index a1cba4fc9..9820b6f9e 100644
--- a/src/algorithms/ctmrg/ctmrg.jl
+++ b/src/algorithms/ctmrg/ctmrg.jl
@@ -88,6 +88,7 @@ containing the following fields:
 
 * `truncation_error` : Last (maximal) SVD truncation error of the CTMRG projectors.
 * `condition_number` : Last (maximal) condition number of the enlarged CTMRG environment.
+* `convergence_error` : Convergence error of the CTMRG algorithm at termination.
 
 In case the `alg` is a `SimultaneousCTMRG`, the last SVD will also be returned:
 
@@ -120,6 +121,7 @@ function leading_boundary(
         for iter in 1:(alg.maxiter)
             env, info = ctmrg_iteration(network, env, alg)  # Grow and renormalize in all 4 directions
             η, CS, TS = calc_convergence(env, CS, TS)
+            info = @insert info.convergence_error = η
 
             if η ≤ alg.tol && iter ≥ alg.miniter
                 ctmrg_logfinish!(log, iter, η, network, env)
diff --git a/test/ctmrg/initialization.jl b/test/ctmrg/initialization.jl
new file mode 100644
index 000000000..111f602d4
--- /dev/null
+++ b/test/ctmrg/initialization.jl
@@ -0,0 +1,51 @@
+using Test
+using TensorKit
+using PEPSKit
+using Random
+
+using MPSKitModels: classical_ising
+
+sd = 12345
+
+# toggle symmetry, but same issue for both
+symmetries = [Z2Irrep, Trivial]
+
+χ = 20
+tol = 1.0e-4
+maxiter = 1000
+verbosity = 2
+trscheme = FixedSpaceTruncation()
+boundary_alg = (;
+    alg = :simultaneous,
+    tol,
+    verbosity,
+    trscheme,
+    maxiter,
+)
+
+@testset "CTMRG environment initialization for critical ising with $S symmetry (#255)" for S in symmetries
+    # initialize
+    T = classical_ising(S)
+    O = T[1]
+    n = InfinitePartitionFunction([O O; O O])
+    Venv = S == Z2Irrep ? Z2Space(0 => χ / 2, 1 => χ / 2) : ℂ^χ
+    P = space(O, 2)
+
+    # random, doesn't converge
+    Random.seed!(sd)
+    env0_rand = initialize_environment(n, RandomInitialization(), Venv)
+    env_rand, info = leading_boundary(env0_rand, n; boundary_alg...)
+    @test info.convergence_error > tol
+
+    # embedded product state, converges
+    Random.seed!(sd)
+    env0_prod = initialize_environment(n, ProductStateInitialization(), Venv)
+    env_prod, info = leading_boundary(env0_prod, n; boundary_alg...)
+    @test info.convergence_error ≤ tol
+
+    # grown product state, converges
+    Random.seed!(sd)
+    env0_appl = initialize_environment(n, ApplicationInitialization(), truncdim(χ))
+    env_appl, info = leading_boundary(env0_appl, n; boundary_alg...)
+    @test info.convergence_error ≤ tol
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 496fb7437..aa153bd73 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -35,6 +35,9 @@ end
         @time @safetestset "correlation length" begin
             include("ctmrg/correlation_length.jl")
         end
+        @time @safetestset "initialization" begin
+            include("ctmrg/initialization.jl")
+        end
     end
     if GROUP == "ALL" || GROUP == "GRADIENTS"
         @time @safetestset "CTMRG gradients" begin

From 2a6f121e6d77796717c37290a438755ae30c7bb5 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Fri, 26 Sep 2025 11:26:22 +0200
Subject: [PATCH 03/24] Not converging isn't really passing though

---
 test/ctmrg/initialization.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/ctmrg/initialization.jl b/test/ctmrg/initialization.jl
index 111f602d4..5ed54fcdb 100644
--- a/test/ctmrg/initialization.jl
+++ b/test/ctmrg/initialization.jl
@@ -35,7 +35,7 @@ boundary_alg = (;
     Random.seed!(sd)
     env0_rand = initialize_environment(n, RandomInitialization(), Venv)
     env_rand, info = leading_boundary(env0_rand, n; boundary_alg...)
-    @test info.convergence_error > tol
+    @test_broken info.convergence_error ≤ tol
 
     # embedded product state, converges
     Random.seed!(sd)

From 78cef7ca319f0c6820fa52021b1b883aab6057aa Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Fri, 26 Sep 2025 13:03:14 +0200
Subject: [PATCH 04/24] Pass virtual space specification through to `CTMRGEnv`
 constructor

---
 src/algorithms/ctmrg/initialization.jl | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/src/algorithms/ctmrg/initialization.jl b/src/algorithms/ctmrg/initialization.jl
index cb5a5eebf..36bd16f15 100644
--- a/src/algorithms/ctmrg/initialization.jl
+++ b/src/algorithms/ctmrg/initialization.jl
@@ -3,23 +3,26 @@ struct ProductStateInitialization <: InitializationStyle end
 struct RandomInitialization <: InitializationStyle end
 struct ApplicationInitialization <: InitializationStyle end
 
+_to_tuple(x) = (x,)
+_to_tuple(x::Tuple) = x
+
 function initialize_environment(
         elt::Type{<:Number},
         n::InfiniteSquareNetwork,
         ::RandomInitialization,
-        init_spec::ElementarySpace = oneunit(spacetype(n)), # TODO: non-uniform space specification?
+        init_spec = oneunit(spacetype(n)),
     )
-    return CTMRGEnv(randn, elt, n, init_spec)
+    return CTMRGEnv(randn, elt, n, _to_tuple(init_spec)...)
 end
 
 function initialize_environment(
         elt::Type{<:Number},
         n::InfiniteSquareNetwork,
         ::ProductStateInitialization,
-        init_spec::ElementarySpace = oneunit(spacetype(n)), # TODO: non-uniform space specification?
+        init_spec = oneunit(spacetype(n)),
     )
     i = one(sectortype(init_spec))
-    env = CTMRGEnv(ones, elt, n, init_spec)
+    env = CTMRGEnv(ones, elt, n, _to_tuple(init_spec)...)
     for (dir, r, c) in Iterators.product(axes(env)...)
         @assert i in blocksectors(env.corners[dir, r, c])
         block(env.corners[dir, r, c], i)[1, 1] = 1

From c20aab3e9d13a7983cc0984ed1369398a6dca529 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Sun, 28 Sep 2025 12:45:49 +0200
Subject: [PATCH 05/24] Slurp, make an actual product state, and increase
 coverage

---
 src/algorithms/ctmrg/initialization.jl | 22 +++++++++++-----------
 test/ctmrg/initialization.jl           |  2 +-
 2 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/algorithms/ctmrg/initialization.jl b/src/algorithms/ctmrg/initialization.jl
index 36bd16f15..14325bfa1 100644
--- a/src/algorithms/ctmrg/initialization.jl
+++ b/src/algorithms/ctmrg/initialization.jl
@@ -3,29 +3,29 @@ struct ProductStateInitialization <: InitializationStyle end
 struct RandomInitialization <: InitializationStyle end
 struct ApplicationInitialization <: InitializationStyle end
 
-_to_tuple(x) = (x,)
-_to_tuple(x::Tuple) = x
-
 function initialize_environment(
         elt::Type{<:Number},
         n::InfiniteSquareNetwork,
         ::RandomInitialization,
-        init_spec = oneunit(spacetype(n)),
+        virtual_spaces... = oneunit(spacetype(n)),
     )
-    return CTMRGEnv(randn, elt, n, _to_tuple(init_spec)...)
+    return CTMRGEnv(randn, elt, n, virtual_spaces...)
 end
 
 function initialize_environment(
         elt::Type{<:Number},
         n::InfiniteSquareNetwork,
         ::ProductStateInitialization,
-        init_spec = oneunit(spacetype(n)),
+        virtual_spaces... = oneunit(spacetype(n)),
     )
-    i = one(sectortype(init_spec))
-    env = CTMRGEnv(ones, elt, n, _to_tuple(init_spec)...)
+    i = one(sectortype(n))
+    env = CTMRGEnv(ones, elt, n, virtual_spaces...)
     for (dir, r, c) in Iterators.product(axes(env)...)
         @assert i in blocksectors(env.corners[dir, r, c])
-        block(env.corners[dir, r, c], i)[1, 1] = 1
+        for (c, b) in blocks(env.corners[dir, r, c])
+            b .= 0
+            c == i && (b[1, 1] = 1)
+        end
     end
     return env
 end
@@ -34,9 +34,9 @@ function initialize_environment(
         elt::Type{<:Number},
         n::InfiniteSquareNetwork,
         ::ApplicationInitialization,
-        init_spec::TruncationScheme;
+        trscheme::TruncationScheme;
         boundary_alg = (;
-            alg = :sequential, tol = 1.0e-5, maxiter = 10, verbosity = -1, trscheme = init_spec,
+            alg = :sequential, tol = 1.0e-5, maxiter = 10, verbosity = -1, trscheme,
         )
     )
     env = initialize_environment(elt, n, ProductStateInitialization())
diff --git a/test/ctmrg/initialization.jl b/test/ctmrg/initialization.jl
index 5ed54fcdb..ba7a2fdd0 100644
--- a/test/ctmrg/initialization.jl
+++ b/test/ctmrg/initialization.jl
@@ -45,7 +45,7 @@ boundary_alg = (;
 
     # grown product state, converges
     Random.seed!(sd)
-    env0_appl = initialize_environment(n, ApplicationInitialization(), truncdim(χ))
+    env0_appl = initialize_environment(InfiniteSquareNetwork(n), ApplicationInitialization(), truncdim(χ))
     env_appl, info = leading_boundary(env0_appl, n; boundary_alg...)
     @test info.convergence_error ≤ tol
 end

From 11d8535e2a34ba1ab57c7fdb1e6ae64e27a25b70 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Sun, 28 Sep 2025 13:19:09 +0200
Subject: [PATCH 06/24] Better optional alg specification

---
 src/algorithms/ctmrg/initialization.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/algorithms/ctmrg/initialization.jl b/src/algorithms/ctmrg/initialization.jl
index 14325bfa1..9308cf128 100644
--- a/src/algorithms/ctmrg/initialization.jl
+++ b/src/algorithms/ctmrg/initialization.jl
@@ -36,9 +36,10 @@ function initialize_environment(
         ::ApplicationInitialization,
         trscheme::TruncationScheme;
         boundary_alg = (;
-            alg = :sequential, tol = 1.0e-5, maxiter = 10, verbosity = -1, trscheme,
+            alg = :sequential, tol = 1.0e-5, maxiter = 10, verbosity = -1,
         )
     )
+    boundary_alg = (; boundary_alg..., trscheme) # merge trscheme with optional alg definition
     env = initialize_environment(elt, n, ProductStateInitialization())
     env, = leading_boundary(env, n; boundary_alg...)
     return env

From 991ecc70abb4f21b56c60a2e238e193f69e4ad45 Mon Sep 17 00:00:00 2001
From: Lander Burgelman <39218680+leburgel@users.noreply.github.com>
Date: Wed, 1 Oct 2025 10:31:07 +0200
Subject: [PATCH 07/24] Apply suggestions from code review

Co-authored-by: Lukas Devos <ldevos98@gmail.com>
---
 src/algorithms/ctmrg/initialization.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/algorithms/ctmrg/initialization.jl b/src/algorithms/ctmrg/initialization.jl
index 9308cf128..0dc0cf65d 100644
--- a/src/algorithms/ctmrg/initialization.jl
+++ b/src/algorithms/ctmrg/initialization.jl
@@ -46,10 +46,10 @@ function initialize_environment(
 end
 
 function initialize_environment(n::InfiniteSquareNetwork, args...; kwargs...)
-    return initialize_environment(ComplexF64, n, args...; kwargs...)
+    return initialize_environment(scalartype(n), n, args...; kwargs...)
 end
 function initialize_environment(A::Union{InfinitePEPS, InfinitePartitionFunction}, args...; kwargs...)
-    return initialize_environment(ComplexF64, A, args...; kwargs...)
+    return initialize_environment(scalartype(A), A, args...; kwargs...)
 end
 function initialize_environment(elt::Type{<:Number}, A::Union{InfinitePEPS, InfinitePartitionFunction}, args...; kwargs...)
     return initialize_environment(elt, InfiniteSquareNetwork(A), args...; kwargs...)

From 4996e1bc56ca025e2a79370124fd4ef8f7f1c486 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Mon, 6 Oct 2025 11:18:28 +0200
Subject: [PATCH 08/24] Add initialization function to
 `ProductStateInitialization` struct

---
 src/algorithms/ctmrg/initialization.jl | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/algorithms/ctmrg/initialization.jl b/src/algorithms/ctmrg/initialization.jl
index 0dc0cf65d..45e2d456e 100644
--- a/src/algorithms/ctmrg/initialization.jl
+++ b/src/algorithms/ctmrg/initialization.jl
@@ -1,15 +1,18 @@
 abstract type InitializationStyle end
 struct ProductStateInitialization <: InitializationStyle end
-struct RandomInitialization <: InitializationStyle end
+struct RandomInitialization{F} <: InitializationStyle
+    f::F
+    RandomInitialization(f::F = randn) where {F} = new{F}(f)
+end
 struct ApplicationInitialization <: InitializationStyle end
 
 function initialize_environment(
         elt::Type{<:Number},
         n::InfiniteSquareNetwork,
-        ::RandomInitialization,
+        alg::RandomInitialization,
         virtual_spaces... = oneunit(spacetype(n)),
     )
-    return CTMRGEnv(randn, elt, n, virtual_spaces...)
+    return CTMRGEnv(alg.f, elt, n, virtual_spaces...)
 end
 
 function initialize_environment(

From f4a9afbe870227160bfa8d001f86b832d6edb62f Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sun, 11 Jan 2026 19:47:02 +0800
Subject: [PATCH 09/24] Change `trscheme` to `trunc`

---
 src/algorithms/ctmrg/initialization.jl | 4 ++--
 test/ctmrg/initialization.jl           | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/algorithms/ctmrg/initialization.jl b/src/algorithms/ctmrg/initialization.jl
index 45e2d456e..c6b3a0d91 100644
--- a/src/algorithms/ctmrg/initialization.jl
+++ b/src/algorithms/ctmrg/initialization.jl
@@ -37,12 +37,12 @@ function initialize_environment(
         elt::Type{<:Number},
         n::InfiniteSquareNetwork,
         ::ApplicationInitialization,
-        trscheme::TruncationScheme;
+        trunc::TruncationStrategy;
         boundary_alg = (;
             alg = :sequential, tol = 1.0e-5, maxiter = 10, verbosity = -1,
         )
     )
-    boundary_alg = (; boundary_alg..., trscheme) # merge trscheme with optional alg definition
+    boundary_alg = (; boundary_alg..., trunc) # merge trunc with optional alg definition
     env = initialize_environment(elt, n, ProductStateInitialization())
     env, = leading_boundary(env, n; boundary_alg...)
     return env
diff --git a/test/ctmrg/initialization.jl b/test/ctmrg/initialization.jl
index ba7a2fdd0..09a99d25b 100644
--- a/test/ctmrg/initialization.jl
+++ b/test/ctmrg/initialization.jl
@@ -14,12 +14,12 @@ symmetries = [Z2Irrep, Trivial]
 tol = 1.0e-4
 maxiter = 1000
 verbosity = 2
-trscheme = FixedSpaceTruncation()
+trunc = FixedSpaceTruncation()
 boundary_alg = (;
     alg = :simultaneous,
     tol,
     verbosity,
-    trscheme,
+    trunc,
     maxiter,
 )
 

From 5a7605bc4034b0a04c1f52e16ff868c460d6a5b3 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Fri, 10 Apr 2026 11:36:07 +0200
Subject: [PATCH 10/24] Stash update

---
 src/PEPSKit.jl                                |   3 +-
 src/algorithms/ctmrg/initialization.jl        |  59 ++++----
 .../product_state_environments.jl             | 141 ++++++++++++++++++
 test/ctmrg/initialization.jl                  |   8 +-
 4 files changed, 175 insertions(+), 36 deletions(-)
 create mode 100644 src/environments/product_state_environments.jl

diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 1a3645c45..e3693cb36 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -68,6 +68,7 @@ include("environments/ctmrg_environments.jl")
 include("environments/vumps_environments.jl")
 include("environments/suweight.jl")
 include("environments/bp_environments.jl")
+include("environments/product_state_environments.jl")
 
 include("algorithms/contractions/ctmrg/types.jl")
 include("algorithms/contractions/ctmrg/expr.jl")
@@ -129,7 +130,7 @@ using .Defaults: set_scheduler!
 export set_scheduler!
 export EighAdjoint, IterEigh, SVDAdjoint, IterSVD, QRAdjoint
 export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG
-export initialize_environment,
+export initialize_ctmrg_environment,
     RandomInitialization, ProductStateInitialization, ApplicationInitialization
 export FixedSpaceTruncation, SiteDependentTruncation
 export HalfInfiniteProjector, FullInfiniteProjector
diff --git a/src/algorithms/ctmrg/initialization.jl b/src/algorithms/ctmrg/initialization.jl
index c6b3a0d91..196e099ac 100644
--- a/src/algorithms/ctmrg/initialization.jl
+++ b/src/algorithms/ctmrg/initialization.jl
@@ -1,12 +1,19 @@
 abstract type InitializationStyle end
-struct ProductStateInitialization <: InitializationStyle end
+struct ProductStateInitialization{F} <: InitializationStyle
+    f::F
+    ProductStateInitialization(f::F = ones) where {F} = new{F}(f)
+end
 struct RandomInitialization{F} <: InitializationStyle
     f::F
     RandomInitialization(f::F = randn) where {F} = new{F}(f)
 end
-struct ApplicationInitialization <: InitializationStyle end
+struct ApplicationInitialization{F} <: InitializationStyle
+    f::F
+    ApplicationInitialization(f::F = ones) where {F} = new{F}(f)
+end
 
-function initialize_environment(
+# initialize randomly, using same virtual space specification as the CTMRGEnv constructor
+function initialize_ctmrg_environment(
         elt::Type{<:Number},
         n::InfiniteSquareNetwork,
         alg::RandomInitialization,
@@ -15,45 +22,35 @@ function initialize_environment(
     return CTMRGEnv(alg.f, elt, n, virtual_spaces...)
 end
 
-function initialize_environment(
+function initialize_ctmrg_environment(
         elt::Type{<:Number},
         n::InfiniteSquareNetwork,
-        ::ProductStateInitialization,
-        virtual_spaces... = oneunit(spacetype(n)),
+        alg::ProductStateInitialization,
     )
-    i = one(sectortype(n))
-    env = CTMRGEnv(ones, elt, n, virtual_spaces...)
-    for (dir, r, c) in Iterators.product(axes(env)...)
-        @assert i in blocksectors(env.corners[dir, r, c])
-        for (c, b) in blocks(env.corners[dir, r, c])
-            b .= 0
-            c == i && (b[1, 1] = 1)
-        end
-    end
+    env = CTMRGEnv(ProductStateEnv(alg.f, elt, n))
     return env
 end
 
-function initialize_environment(
+function initialize_ctmrg_environment(
         elt::Type{<:Number},
         n::InfiniteSquareNetwork,
-        ::ApplicationInitialization,
-        trunc::TruncationStrategy;
-        boundary_alg = (;
-            alg = :sequential, tol = 1.0e-5, maxiter = 10, verbosity = -1,
-        )
+        alg::ApplicationInitialization,
+        env0::ProductStateEnv = ProductStateEnv(alg.f, elt, n)
     )
-    boundary_alg = (; boundary_alg..., trunc) # merge trunc with optional alg definition
-    env = initialize_environment(elt, n, ProductStateInitialization())
-    env, = leading_boundary(env, n; boundary_alg...)
+    dummy_alg = CTMRGAlgorithm(; alg = :simultaneous, trunc = (; alg = :notrunc))
+    env, = ctmrg_iteration(n, CTMRGEnv(env0), dummy_alg)
     return env
 end
 
-function initialize_environment(n::InfiniteSquareNetwork, args...; kwargs...)
-    return initialize_environment(scalartype(n), n, args...; kwargs...)
-end
-function initialize_environment(A::Union{InfinitePEPS, InfinitePartitionFunction}, args...; kwargs...)
-    return initialize_environment(scalartype(A), A, args...; kwargs...)
+function initialize_ctmrg_environment(
+        A::Union{InfiniteSquareNetwork, InfinitePEPS, InfinitePartitionFunction}, args...;
+        kwargs...
+    )
+    return initialize_ctmrg_environment(scalartype(A), A, args...; kwargs...)
 end
-function initialize_environment(elt::Type{<:Number}, A::Union{InfinitePEPS, InfinitePartitionFunction}, args...; kwargs...)
-    return initialize_environment(elt, InfiniteSquareNetwork(A), args...; kwargs...)
+function initialize_ctmrg_environment(
+        elt::Type{<:Number}, A::Union{InfinitePEPS, InfinitePartitionFunction}, args...;
+        kwargs...
+    )
+    return initialize_ctmrg_environment(elt, InfiniteSquareNetwork(A), args...; kwargs...)
 end
diff --git a/src/environments/product_state_environments.jl b/src/environments/product_state_environments.jl
new file mode 100644
index 000000000..7e009c5c4
--- /dev/null
+++ b/src/environments/product_state_environments.jl
@@ -0,0 +1,141 @@
+"""
+$(TYPEDEF)
+
+Tensor product environment for an infinite square network, containing a 4 x rows x cols
+array of tensors, defined for each nearest neighbor bond in the network. 
+
+The product state tensors `p` connect to the network tensors 
+`P` at site `[r,c]` in the unit cell as:
+```
+                    p[1,r-1,c]
+                    |
+    p[4,r,c-1]------P[r,c]------p[2,r,c+1]
+                    |
+                    p[3,r+1,c]
+```
+## Fields
+
+$(TYPEDFIELDS)
+"""
+struct ProductStateEnv{T}
+    "4 x rows x cols array of edge tensors making up a product state environment, where the
+    first dimension specifies the spatial direction"
+    edges::Array{T, 3}
+end
+
+"""
+    ProductStateEnv(
+        [f=randn, T=ComplexF64], Ds_north::A, Ds_east::A
+    ) where {A <: AbstractMatrix{<:ProductSpace}}
+
+Construct a product state environment by specifying matrices of north and east virtual spaces of the
+corresponding [`InfiniteSquareNetwork`](@ref). Each matrix entry corresponds to a site in the unit cell.
+
+Each entry of the `Ds_north` and `Ds_east` matrices corresponds to an effective local space
+of the network, and can be represented as a `ProductSpace` (e.g.
+for the case of a network representing overlaps of PEPSs).
+"""
+function ProductStateEnv(
+        Ds_north::A, Ds_east::A
+    ) where {A <: AbstractMatrix{<:ProductSpace}}
+    return ProductStateEnv(randn, ComplexF64, N, Ds_north, Ds_east)
+end
+function ProductStateEnv(
+        f, T, Ds_north::A, Ds_east::A
+    ) where {A <: AbstractMatrix{<:ProductSpace}}
+    Ds_south = _elementwise_dual.(circshift(Ds_north, (-1, 0)))
+    Ds_west = _elementwise_dual.(circshift(Ds_east, (0, 1)))
+    edges = map(Iterators.product(1:4, axes(Ds_north, 1), axes(Ds_north, 2))) do (dir, r, c)
+        msg = if dir == NORTH
+            f(T, Ds_north[_next(r, end), c])
+        elseif dir == EAST
+            f(T, Ds_east[r, _prev(c, end)])
+        elseif dir == SOUTH
+            f(T, Ds_south[_prev(r, end), c])
+        else # WEST
+            f(T, Ds_west[r, _next(c, end)])
+        end
+        return msg
+    end
+    normalize!.(edges)
+    return ProductStateEnv(edges)
+end
+
+"""
+    ProductStateEnv(
+        [f=randn, T=ComplexF64], D_north::P, D_east::P;
+        unitcell::Tuple{Int, Int} = (1, 1)
+    ) where {P <: ProductSpace}
+
+Construct a product state environment by specifying the north and east virtual spaces of the
+corresponding [`InfiniteSquareNetwork`](@ref). The network unit cell can be specified
+by the `unitcell` keyword argument.
+"""
+function ProductStateEnv(
+        D_north::P, D_east::P;
+        unitcell::Tuple{Int, Int} = (1, 1)
+    ) where {P <: ProductSpace}
+    return ProductStateEnv(randn, ComplexF64, D_north, D_east; unitcell)
+end
+function ProductStateEnv(
+        f, T, D_north::P, D_east::P;
+        unitcell::Tuple{Int, Int} = (1, 1)
+    ) where {P <: ProductSpace}
+    return ProductStateEnv(f, T, N, fill(D_north, unitcell), fill(D_east, unitcell))
+end
+
+"""
+    ProductStateEnv([f=ones, T=ComplexF64], network::InfiniteSquareNetwork)
+
+Construct a product state environment by specifying a corresponding [`InfiniteSquareNetwork`](@ref).
+"""
+function ProductStateEnv(f, T, network::InfiniteSquareNetwork)
+    Ds_north = _north_edge_physical_spaces(network)
+    Ds_east = _east_edge_physical_spaces(network)
+    return ProductStateEnv(f, T, Ds_north, Ds_east)
+end
+function ProductStateEnv(network::InfiniteSquareNetwork)
+    return ProductStateEnv(ones, scalartype(network), network) # TODO: do we want to use a different default function?
+end
+
+function ProductStateEnv(state::Union{InfinitePartitionFunction, InfinitePEPS, InfinitePEPO}, args...; kwargs...)
+    return ProductStateEnv(InfiniteSquareNetwork(state), args...; kwargs...)
+end
+function ProductStateEnv(state::Union{InfinitePEPS, InfinitePEPO}, args...; kwargs...)
+    return ProductStateEnv(InfiniteSquareNetwork(state), args...; kwargs...)
+end
+function ProductStateEnv(f, T, state::Union{InfinitePartitionFunction, InfinitePEPS, InfinitePEPO}, args...; kwargs...)
+    return ProductStateEnv(f, T, InfiniteSquareNetwork(state), args...; kwargs...)
+end
+
+Base.eltype(::Type{ProductStateEnv{T}}) where {T} = T
+Base.size(env::ProductStateEnv, args...) = size(env.edges, args...)
+Base.getindex(env::ProductStateEnv, args...) = Base.getindex(env.edges, args...)
+Base.axes(env::ProductStateEnv, args...) = Base.axes(env.edges, args...)
+Base.eachindex(env::ProductStateEnv) = eachindex(IndexCartesian(), env.edges)
+VectorInterface.scalartype(::Type{ProductStateEnv{T}}) where {T} = scalartype(T)
+TensorKit.spacetype(::Type{ProductStateEnv{T}}) where {T} = spacetype(T)
+
+function eachcoordinate(x::ProductStateEnv)
+    return collect(Iterators.product(axes(x, 2), axes(x, 3)))
+end
+function eachcoordinate(x::ProductStateEnv, dirs)
+    return collect(Iterators.product(dirs, axes(x, 2), axes(x, 3)))
+end
+
+# conversion to CTMRGEnv
+"""
+    CTMRGEnv(prod_env::ProductStateEnv)
+
+Construct a CTMRG environment with a trivial virtual space of bond dimension χ = 1 
+from the product state environment `prod_env`.
+"""
+function CTMRGEnv(prod_env::ProductStateEnv)
+    edges = map(CartesianIndices(prod_env.edges)) do idx
+        return insertleftunit(insertleftunit(prod_env.edges[idx]), 1)
+    end
+    corners = map(CartesianIndices(edges)) do _
+        return TensorKit.id(scalartype(prod_env), oneunit(spacetype(prod_env)))
+    end
+    return CTMRGEnv(corners, edges)
+end
diff --git a/test/ctmrg/initialization.jl b/test/ctmrg/initialization.jl
index 09a99d25b..14b77bf67 100644
--- a/test/ctmrg/initialization.jl
+++ b/test/ctmrg/initialization.jl
@@ -14,7 +14,7 @@ symmetries = [Z2Irrep, Trivial]
 tol = 1.0e-4
 maxiter = 1000
 verbosity = 2
-trunc = FixedSpaceTruncation()
+trunc = truncrank(χ)
 boundary_alg = (;
     alg = :simultaneous,
     tol,
@@ -33,19 +33,19 @@ boundary_alg = (;
 
     # random, doesn't converge
     Random.seed!(sd)
-    env0_rand = initialize_environment(n, RandomInitialization(), Venv)
+    env0_rand = initialize_ctmrg_environment(n, RandomInitialization(), Venv)
     env_rand, info = leading_boundary(env0_rand, n; boundary_alg...)
     @test_broken info.convergence_error ≤ tol
 
     # embedded product state, converges
     Random.seed!(sd)
-    env0_prod = initialize_environment(n, ProductStateInitialization(), Venv)
+    env0_prod = initialize_ctmrg_environment(n, ProductStateInitialization())
     env_prod, info = leading_boundary(env0_prod, n; boundary_alg...)
     @test info.convergence_error ≤ tol
 
     # grown product state, converges
     Random.seed!(sd)
-    env0_appl = initialize_environment(InfiniteSquareNetwork(n), ApplicationInitialization(), truncdim(χ))
+    env0_appl = initialize_ctmrg_environment(n, ApplicationInitialization())
     env_appl, info = leading_boundary(env0_appl, n; boundary_alg...)
     @test info.convergence_error ≤ tol
 end

From 25026152b39ffa6beea2ecb3f59b1d73b38ca036 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Sat, 11 Apr 2026 10:52:20 +0800
Subject: [PATCH 11/24] Update env init in finite-T SU tests

---
 test/timeevol/j1j2_finiteT.jl     | 6 +++---
 test/timeevol/tf_ising_finiteT.jl | 8 +++-----
 2 files changed, 6 insertions(+), 8 deletions(-)

diff --git a/test/timeevol/j1j2_finiteT.jl b/test/timeevol/j1j2_finiteT.jl
index 63576e6a1..c22c05251 100644
--- a/test/timeevol/j1j2_finiteT.jl
+++ b/test/timeevol/j1j2_finiteT.jl
@@ -10,9 +10,9 @@ using PEPSKit
 bm = [-0.1235, -0.213]
 
 function converge_env(state, χ::Int)
-    trunc1 = truncrank(χ) & truncerror(; atol = 1.0e-12)
-    env0 = CTMRGEnv(ones, Float64, state, Vect[SU2Irrep](0 => 1))
-    env, = leading_boundary(env0, state; alg = :sequential, trunc = trunc1, tol = 1.0e-10)
+    env0 = initialize_ctmrg_environment(state, ProductStateInitialization())
+    trunc = truncrank(χ) & truncerror(; atol = 1.0e-12)
+    env, = leading_boundary(env0, state; alg = :sequential, trunc, tol = 1.0e-10)
     return env
 end
 
diff --git a/test/timeevol/tf_ising_finiteT.jl b/test/timeevol/tf_ising_finiteT.jl
index 6164460f1..fa664da9e 100644
--- a/test/timeevol/tf_ising_finiteT.jl
+++ b/test/timeevol/tf_ising_finiteT.jl
@@ -10,11 +10,9 @@ bm_β = [0.5632, 0.0]
 bm_2β = [0.5297, 0.8265]
 
 function converge_env(state, χ::Int)
-    trunc1 = truncrank(4) & truncerror(; atol = 1.0e-12)
-    env0 = CTMRGEnv(rand, Float64, state, ℂ^2)
-    env, = leading_boundary(env0, state; alg = :sequential, trunc = trunc1, tol = 1.0e-10)
-    trunc2 = truncrank(χ) & truncerror(; atol = 1.0e-12)
-    env, = leading_boundary(env, state; alg = :sequential, trunc = trunc2, tol = 1.0e-10)
+    env0 = initialize_ctmrg_environment(state, ProductStateInitialization())
+    trunc = truncrank(χ) & truncerror(; atol = 1.0e-12)
+    env, = leading_boundary(env0, state; alg = :sequential, trunc, tol = 1.0e-10)
     return env
 end
 

From 885faaf16a9b2938df7438b7d345832f8fe3c534 Mon Sep 17 00:00:00 2001
From: Lander Burgelman <39218680+leburgel@users.noreply.github.com>
Date: Wed, 15 Apr 2026 10:52:00 +0200
Subject: [PATCH 12/24] Update test/ctmrg/initialization.jl

Co-authored-by: Yue Zhengyuan <yuezy1997@icloud.com>
---
 test/ctmrg/initialization.jl | 8 +-------
 1 file changed, 1 insertion(+), 7 deletions(-)

diff --git a/test/ctmrg/initialization.jl b/test/ctmrg/initialization.jl
index 14b77bf67..ea8a4ed25 100644
--- a/test/ctmrg/initialization.jl
+++ b/test/ctmrg/initialization.jl
@@ -15,13 +15,7 @@ tol = 1.0e-4
 maxiter = 1000
 verbosity = 2
 trunc = truncrank(χ)
-boundary_alg = (;
-    alg = :simultaneous,
-    tol,
-    verbosity,
-    trunc,
-    maxiter,
-)
+boundary_alg = (; alg = :simultaneous, tol, verbosity, trunc, maxiter)
 
 @testset "CTMRG environment initialization for critical ising with $S symmetry (#255)" for S in symmetries
     # initialize

From d56074aabeb01ec0a83ff2e28d9cec83472ca63a Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Fri, 5 Jun 2026 13:32:23 +0200
Subject: [PATCH 13/24] Some fixes, and add environment initialization for PEPS
 using identity edges

---
 src/PEPSKit.jl                         |  2 +-
 src/algorithms/ctmrg/ctmrg.jl          |  1 -
 src/algorithms/ctmrg/initialization.jl |  2 +-
 src/utility/util.jl                    | 14 ++++++++
 test/ctmrg/initialization.jl           | 47 ++++++++++++++++++++++++--
 5 files changed, 60 insertions(+), 6 deletions(-)

diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 7f4edd07c..7f534d6f9 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -144,7 +144,7 @@ export HalfInfiniteProjector, FullInfiniteProjector
 export C4vCTMRG, C4vEighProjector, C4vQRProjector
 export initialize_random_c4v_env, initialize_singlet_c4v_env
 export LocalOperator, physicalspace
-export product_peps
+export product_peps, bipartite_id
 export reduced_densitymatrix, expectation_value, network_value, cost_function
 export correlator, correlation_length
 export leading_boundary
diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl
index 4e2bfc2ce..57e7a901d 100644
--- a/src/algorithms/ctmrg/ctmrg.jl
+++ b/src/algorithms/ctmrg/ctmrg.jl
@@ -124,7 +124,6 @@ function leading_boundary(
         for iter in 1:(alg.maxiter)
             env, info_iter = ctmrg_iteration(network, env, alg)
             η, CS, TS = calc_convergence(env, CS, TS)
-            info = @insert info.convergence_error = η
 
             if η ≤ alg.tol && iter ≥ alg.miniter
                 ctmrg_logfinish!(log, iter, η, network, env)
diff --git a/src/algorithms/ctmrg/initialization.jl b/src/algorithms/ctmrg/initialization.jl
index 196e099ac..d8dedf094 100644
--- a/src/algorithms/ctmrg/initialization.jl
+++ b/src/algorithms/ctmrg/initialization.jl
@@ -37,7 +37,7 @@ function initialize_ctmrg_environment(
         alg::ApplicationInitialization,
         env0::ProductStateEnv = ProductStateEnv(alg.f, elt, n)
     )
-    dummy_alg = CTMRGAlgorithm(; alg = :simultaneous, trunc = (; alg = :notrunc))
+    dummy_alg = SimultaneousCTMRG(trunc = (; alg = :notrunc))
     env, = ctmrg_iteration(n, CTMRGEnv(env0), dummy_alg)
     return env
 end
diff --git a/src/utility/util.jl b/src/utility/util.jl
index 2c58f8801..f6d071bc8 100644
--- a/src/utility/util.jl
+++ b/src/utility/util.jl
@@ -208,3 +208,17 @@ function _permute_to_last(axes::NTuple{N, Int}, ax::Int) where {N}
     new_axes = (ntuple(i -> axes[biperm[1][i]], N - 1)..., ax)
     return new_axes, biperm
 end
+
+"""
+    bipartite_id([T::Type=Float64], V::ProductSpace{S, 2}) where {S}
+
+Constructs a tensor corresonding to the (permutation of) the identity operator between
+two componens of a bipartite product space.
+"""
+bipartite_id(V::ProductSpace) = bipartite_id(Float64, V)
+function bipartite_id(T::Type, V::ProductSpace)
+    throw(ArgumentError("bipartite_id is only defined for ProductSpace with 2 components, but got $(length(V)) components."))
+end
+function bipartite_id(T::Type, V::ProductSpace{S, 2}) where {S}
+    return permute(isomorphism(T, V[1] ← V[2]'), ((1, 2), ()))
+end
diff --git a/test/ctmrg/initialization.jl b/test/ctmrg/initialization.jl
index ea8a4ed25..29672d8bb 100644
--- a/test/ctmrg/initialization.jl
+++ b/test/ctmrg/initialization.jl
@@ -9,20 +9,24 @@ sd = 12345
 
 # toggle symmetry, but same issue for both
 symmetries = [Z2Irrep, Trivial]
+make_space(::Type{Z2Irrep}, d::Int) = Z2Space(0 => d / 2, 1 => d / 2)
+make_space(::Type{Trivial}, d::Int) = ComplexSpace(d)
 
+d = 2
+D = 4
 χ = 20
 tol = 1.0e-4
 maxiter = 1000
 verbosity = 2
 trunc = truncrank(χ)
-boundary_alg = (; alg = :simultaneous, tol, verbosity, trunc, maxiter)
+boundary_alg = (; alg = :SimultaneousCTMRG, tol, verbosity, trunc, maxiter)
 
 @testset "CTMRG environment initialization for critical ising with $S symmetry (#255)" for S in symmetries
     # initialize
     T = classical_ising(S)
     O = T[1]
     n = InfinitePartitionFunction([O O; O O])
-    Venv = S == Z2Irrep ? Z2Space(0 => χ / 2, 1 => χ / 2) : ℂ^χ
+    Venv = make_space(S, χ)
     P = space(O, 2)
 
     # random, doesn't converge
@@ -31,7 +35,7 @@ boundary_alg = (; alg = :simultaneous, tol, verbosity, trunc, maxiter)
     env_rand, info = leading_boundary(env0_rand, n; boundary_alg...)
     @test_broken info.convergence_error ≤ tol
 
-    # embedded product state, converges
+    # embedded random product state, converges
     Random.seed!(sd)
     env0_prod = initialize_ctmrg_environment(n, ProductStateInitialization())
     env_prod, info = leading_boundary(env0_prod, n; boundary_alg...)
@@ -42,4 +46,41 @@ boundary_alg = (; alg = :simultaneous, tol, verbosity, trunc, maxiter)
     env0_appl = initialize_ctmrg_environment(n, ApplicationInitialization())
     env_appl, info = leading_boundary(env0_appl, n; boundary_alg...)
     @test info.convergence_error ≤ tol
+
+    # PEPS-specific identity initialization; should throw when used on partition functions
+    Random.seed!(sd)
+    @test_throws ArgumentError env0_prod_id = initialize_ctmrg_environment(n, ProductStateInitialization(bipartite_id))
+end
+
+@testset "CTMRG environment initialization for PEPS with $S symmetry (#255)" for S in symmetries
+    # initialize
+    P = make_space(S, d)
+    Vpeps = make_space(S, D)
+    Venv = make_space(S, χ)
+    peps = InfinitePEPS(P, Vpeps; unitcell = (2, 2))
+    n = InfiniteSquareNetwork(peps)
+
+    # random, converges
+    Random.seed!(sd)
+    env0_rand = initialize_ctmrg_environment(n, RandomInitialization(), Venv)
+    env_rand, info = leading_boundary(env0_rand, n; boundary_alg...)
+    @test info.convergence_error ≤ tol
+
+    # embedded random product state, converges
+    Random.seed!(sd)
+    env0_prod = initialize_ctmrg_environment(n, ProductStateInitialization())
+    env_prod, info = leading_boundary(env0_prod, n; boundary_alg...)
+    @test info.convergence_error ≤ tol
+
+    # embedded product state as identity from ket to bra, converges
+    Random.seed!(sd)
+    env0_prod = initialize_ctmrg_environment(n, ProductStateInitialization(bipartite_id))
+    env_prod, info = leading_boundary(env0_prod, n; boundary_alg...)
+    @test info.convergence_error ≤ tol
+
+    # grown product state, converges
+    Random.seed!(sd)
+    env0_appl = initialize_ctmrg_environment(n, ApplicationInitialization())
+    env_appl, info = leading_boundary(env0_appl, n; boundary_alg...)
+    @test info.convergence_error ≤ tol
 end

From fa0d980fe434a34b6e18f373cd4196b5a441059c Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Fri, 5 Jun 2026 13:34:50 +0200
Subject: [PATCH 14/24] No longer use `@insert`

---
 src/PEPSKit.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 7f534d6f9..d55c13f12 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -3,7 +3,7 @@ module PEPSKit
 using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf
 using Random
 using Compat
-using Accessors: @set, @reset, @insert
+using Accessors: @set, @reset
 using VectorInterface
 import VectorInterface as VI
 

From 203266a04dd7b9be2f8139028caa92874cf213e1 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@ugent.be>
Date: Mon, 15 Jun 2026 11:53:20 +0200
Subject: [PATCH 15/24] Remove `bipartite_id`, add dedicated initialization
 from identity, add docstrings

---
 src/PEPSKit.jl                         |  5 +-
 src/algorithms/ctmrg/initialization.jl | 69 ++++++++++++++++++++------
 src/utility/initialization.jl          | 67 +++++++++++++++++++++++++
 src/utility/util.jl                    | 14 ------
 test/ctmrg/initialization.jl           |  8 +--
 5 files changed, 127 insertions(+), 36 deletions(-)
 create mode 100644 src/utility/initialization.jl

diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index d55c13f12..a6e199c54 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -50,6 +50,7 @@ include("utility/rotations.jl")
 include("utility/hook_pullback.jl")
 include("utility/autoopt.jl")
 include("utility/retractions.jl")
+include("utility/initialization.jl")
 
 include("networks/tensors.jl")
 include("networks/local_sandwich.jl")
@@ -137,14 +138,14 @@ export set_scheduler!
 export EighAdjoint, IterEigh, SVDAdjoint, IterSVD, QRAdjoint
 export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG
 export initialize_ctmrg_environment,
-    RandomInitialization, ProductStateInitialization, ApplicationInitialization
+    RandomInitialization, ProductStateInitialization, ApplicationInitialization, IdentityInitialization
 export corner, edge, setcorner!, setedge!
 export FixedSpaceTruncation, SiteDependentTruncation
 export HalfInfiniteProjector, FullInfiniteProjector
 export C4vCTMRG, C4vEighProjector, C4vQRProjector
 export initialize_random_c4v_env, initialize_singlet_c4v_env
 export LocalOperator, physicalspace
-export product_peps, bipartite_id
+export product_peps
 export reduced_densitymatrix, expectation_value, network_value, cost_function
 export correlator, correlation_length
 export leading_boundary
diff --git a/src/algorithms/ctmrg/initialization.jl b/src/algorithms/ctmrg/initialization.jl
index d8dedf094..73e7b0769 100644
--- a/src/algorithms/ctmrg/initialization.jl
+++ b/src/algorithms/ctmrg/initialization.jl
@@ -1,18 +1,9 @@
-abstract type InitializationStyle end
-struct ProductStateInitialization{F} <: InitializationStyle
-    f::F
-    ProductStateInitialization(f::F = ones) where {F} = new{F}(f)
-end
-struct RandomInitialization{F} <: InitializationStyle
-    f::F
-    RandomInitialization(f::F = randn) where {F} = new{F}(f)
-end
-struct ApplicationInitialization{F} <: InitializationStyle
-    f::F
-    ApplicationInitialization(f::F = ones) where {F} = new{F}(f)
-end
+"""
+    initialize_ctmrg_environment([elt::Type{<:Number},] n::InfiniteSquareNetwork, alg::RandomInitialization, virtual_spaces...)
 
-# initialize randomly, using same virtual space specification as the CTMRGEnv constructor
+Initialize a fully random `CTMRGEnv` using the given environment virtual spaces. See
+[`CTMRGEnv`](@ref) for details on the expected format of the virtual spaces.
+"""
 function initialize_ctmrg_environment(
         elt::Type{<:Number},
         n::InfiniteSquareNetwork,
@@ -22,6 +13,12 @@ function initialize_ctmrg_environment(
     return CTMRGEnv(alg.f, elt, n, virtual_spaces...)
 end
 
+"""
+    initialize_ctmrg_environment([elt::Type{<:Number},] n::InfiniteSquareNetwork, alg::RandomInitialization)
+
+Initialize a `CTMRGEnv` corresponding to a product state with trivial virtual spaces and
+corners. The product state edge tensors are initialized as `alg.f(elt, V::ProductSpace)`.
+"""
 function initialize_ctmrg_environment(
         elt::Type{<:Number},
         n::InfiniteSquareNetwork,
@@ -31,14 +28,54 @@ function initialize_ctmrg_environment(
     return env
 end
 
+_CTMRGEnv(env) = CTMRGEnv(env)
+_CTMRGEnv(env::CTMRGEnv) = env
+
+"""
+    initialize_ctmrg_environment([elt::Type{<:Number},] n::InfiniteSquareNetwork, alg::RandomInitialization, [env0])
+
+Initialize a `CTMRGEnv` by applying a single untruncated iteration of
+[`SimultaneousCTMRG`](@ref) to a given initial environment. By default, the starting
+environment is chosen as a random product state.
+"""
 function initialize_ctmrg_environment(
         elt::Type{<:Number},
         n::InfiniteSquareNetwork,
         alg::ApplicationInitialization,
-        env0::ProductStateEnv = ProductStateEnv(alg.f, elt, n)
+        env0 = ProductStateEnv(alg.f, elt, n)
     )
     dummy_alg = SimultaneousCTMRG(trunc = (; alg = :notrunc))
-    env, = ctmrg_iteration(n, CTMRGEnv(env0), dummy_alg)
+    env, = ctmrg_iteration(n, _CTMRGEnv(env0), dummy_alg)
+    return env
+end
+
+_check_two_layer(::InfiniteSquareNetwork) = false
+_check_two_layer(::InfiniteSquareNetwork{<:PEPSSandwich}) = true
+
+"""
+    initialize_ctmrg_environment([elt::Type{<:Number},] n::InfiniteSquareNetwork, alg::RandomInitialization, [env0])
+
+Initialize a `CTMRGEnv` corresponding to a product state acting as an identity between the
+virtual spaces of a two-layer network, for example
+```
+         ╱      
+┌-----ket----- 
+|    ╱ |        
+|      |   
+|      | ╱      
+└-----bra----- 
+     ╱          
+```
+"""
+function initialize_ctmrg_environment(
+        elt::Type{<:Number},
+        n::InfiniteSquareNetwork,
+        ::IdentityInitialization,
+    )
+    _check_two_layer(n) ||
+        throw(ArgumentError("Identity initialization is only defined for two-layer networks."))
+    bp_env = BPEnv(isomorphism, elt, n)
+    env = CTMRGEnv(bp_env)
     return env
 end
 
diff --git a/src/utility/initialization.jl b/src/utility/initialization.jl
new file mode 100644
index 000000000..cf44ef463
--- /dev/null
+++ b/src/utility/initialization.jl
@@ -0,0 +1,67 @@
+"""
+$(TYPEDEF)
+
+Abstract super type for different initialization strategies for contraction environments.
+"""
+abstract type InitializationStyle end
+
+"""
+$(TYPEDEF)
+
+Initialize a contraction environment from a product state made up of `(N, 0)` tensors.
+
+## Constructors
+
+    ProductStateInitialization(f = ones)
+
+Contructs a product state initialization strategy, where the product state tensors
+are initialized by the function `f` as `f(T::Type{<:Number}, V::ProductSpace)`.
+"""
+struct ProductStateInitialization{F} <: InitializationStyle
+    f::F
+    ProductStateInitialization(f::F = ones) where {F} = new{F}(f)
+end
+
+"""
+$(TYPEDEF)
+
+Initialize a fully random contraction environment.
+
+## Constructors
+
+    RandomInitialization(f = randn)
+
+Contructs a random initialization strategy, where the environment tensors are initialized by
+the function `f` as `f(T::Type{<:Number}, V::HomSpace)`.
+"""
+struct RandomInitialization{F} <: InitializationStyle
+    f::F
+    RandomInitialization(f::F = randn) where {F} = new{F}(f)
+end
+
+"""
+$(TYPEDEF)
+
+Initialize a contraction environment by applying a single iteration of a contraction
+algorithm to a given environment.
+
+## Constructors
+
+    ApplicationInitialization(f = ones)
+
+Contructs an application initialization strategy, where by default the starting environment
+is initialized using a `ProductStateInitialization(f)` strategy.
+"""
+struct ApplicationInitialization{F} <: InitializationStyle
+    f::F
+    ApplicationInitialization(f::F = ones) where {F} = new{F}(f)
+end
+
+"""
+$(TYPEDEF)
+
+Initialize a contraction environment 
+
+Only works in very specific cases.
+"""
+struct IdentityInitialization <: InitializationStyle end
diff --git a/src/utility/util.jl b/src/utility/util.jl
index f6d071bc8..2c58f8801 100644
--- a/src/utility/util.jl
+++ b/src/utility/util.jl
@@ -208,17 +208,3 @@ function _permute_to_last(axes::NTuple{N, Int}, ax::Int) where {N}
     new_axes = (ntuple(i -> axes[biperm[1][i]], N - 1)..., ax)
     return new_axes, biperm
 end
-
-"""
-    bipartite_id([T::Type=Float64], V::ProductSpace{S, 2}) where {S}
-
-Constructs a tensor corresonding to the (permutation of) the identity operator between
-two componens of a bipartite product space.
-"""
-bipartite_id(V::ProductSpace) = bipartite_id(Float64, V)
-function bipartite_id(T::Type, V::ProductSpace)
-    throw(ArgumentError("bipartite_id is only defined for ProductSpace with 2 components, but got $(length(V)) components."))
-end
-function bipartite_id(T::Type, V::ProductSpace{S, 2}) where {S}
-    return permute(isomorphism(T, V[1] ← V[2]'), ((1, 2), ()))
-end
diff --git a/test/ctmrg/initialization.jl b/test/ctmrg/initialization.jl
index 29672d8bb..0969cbf46 100644
--- a/test/ctmrg/initialization.jl
+++ b/test/ctmrg/initialization.jl
@@ -49,10 +49,10 @@ boundary_alg = (; alg = :SimultaneousCTMRG, tol, verbosity, trunc, maxiter)
 
     # PEPS-specific identity initialization; should throw when used on partition functions
     Random.seed!(sd)
-    @test_throws ArgumentError env0_prod_id = initialize_ctmrg_environment(n, ProductStateInitialization(bipartite_id))
+    @test_throws ArgumentError env0_prod_id = initialize_ctmrg_environment(n, IdentityInitialization())
 end
 
-@testset "CTMRG environment initialization for PEPS with $S symmetry (#255)" for S in symmetries
+@testset "CTMRG environment initialization for PEPS with $S symmetry" for S in symmetries
     # initialize
     P = make_space(S, d)
     Vpeps = make_space(S, D)
@@ -74,8 +74,8 @@ end
 
     # embedded product state as identity from ket to bra, converges
     Random.seed!(sd)
-    env0_prod = initialize_ctmrg_environment(n, ProductStateInitialization(bipartite_id))
-    env_prod, info = leading_boundary(env0_prod, n; boundary_alg...)
+    env0_prod_id = initialize_ctmrg_environment(n, IdentityInitialization())
+    env_prod, info = leading_boundary(env0_prod_id, n; boundary_alg...)
     @test info.convergence_error ≤ tol
 
     # grown product state, converges

From 4c844f9e7cd55b8eaa8c51f476a816cbd46c3c47 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@ugent.be>
Date: Tue, 16 Jun 2026 14:33:20 +0200
Subject: [PATCH 16/24] Reduce boilerplate a bit, add test for starting from
 specific product state

---
 .../product_state_environments.jl             | 77 +++++++------------
 test/ctmrg/initialization.jl                  |  9 +++
 2 files changed, 35 insertions(+), 51 deletions(-)

diff --git a/src/environments/product_state_environments.jl b/src/environments/product_state_environments.jl
index 7e009c5c4..fd660d659 100644
--- a/src/environments/product_state_environments.jl
+++ b/src/environments/product_state_environments.jl
@@ -21,6 +21,21 @@ struct ProductStateEnv{T}
     "4 x rows x cols array of edge tensors making up a product state environment, where the
     first dimension specifies the spatial direction"
     edges::Array{T, 3}
+    ProductStateEnv{T}(edges::Array{T, 3}) where {T} = new{T}(edges)
+    function ProductStateEnv(edges::Array{T, 3}) where {T}
+        foreach(Iterators.product(axes(edges)[2:3]...)) do (d, w)
+            codomain(edges[NORTH, d, w]) == _elementwise_dual(codomain(edges[SOUTH, _prev(d, end), w])) ||
+                throw(
+                SpaceMismatch("North virtual space at site $((d, w)) does not match: $(space(edges[NORTH, d, w])) vs $(space(edges[SOUTH, _prev(d, end), w])).")
+            )
+            codomain(edges[EAST, d, w]) == _elementwise_dual(codomain(edges[WEST, d, _next(w, end)])) ||
+                throw(SpaceMismatch("East virtual space at site $((d, w)) does not match: $(space(edges[EAST, d, w])) vs $(space(edges[WEST, d, _next(w, end)]))."))
+        end
+        foreach(Iterators.product(axes(edges)...)) do (dir, d, w)
+            dim(space(edges[dir, d, w])) > 0 || @warn "no fusion channels for edge ($dir, $d, $w)"
+        end
+        return new{T}(edges)
+    end
 end
 
 """
@@ -35,11 +50,6 @@ Each entry of the `Ds_north` and `Ds_east` matrices corresponds to an effective
 of the network, and can be represented as a `ProductSpace` (e.g.
 for the case of a network representing overlaps of PEPSs).
 """
-function ProductStateEnv(
-        Ds_north::A, Ds_east::A
-    ) where {A <: AbstractMatrix{<:ProductSpace}}
-    return ProductStateEnv(randn, ComplexF64, N, Ds_north, Ds_east)
-end
 function ProductStateEnv(
         f, T, Ds_north::A, Ds_east::A
     ) where {A <: AbstractMatrix{<:ProductSpace}}
@@ -60,32 +70,12 @@ function ProductStateEnv(
     normalize!.(edges)
     return ProductStateEnv(edges)
 end
-
-"""
-    ProductStateEnv(
-        [f=randn, T=ComplexF64], D_north::P, D_east::P;
-        unitcell::Tuple{Int, Int} = (1, 1)
-    ) where {P <: ProductSpace}
-
-Construct a product state environment by specifying the north and east virtual spaces of the
-corresponding [`InfiniteSquareNetwork`](@ref). The network unit cell can be specified
-by the `unitcell` keyword argument.
-"""
-function ProductStateEnv(
-        D_north::P, D_east::P;
-        unitcell::Tuple{Int, Int} = (1, 1)
-    ) where {P <: ProductSpace}
-    return ProductStateEnv(randn, ComplexF64, D_north, D_east; unitcell)
-end
-function ProductStateEnv(
-        f, T, D_north::P, D_east::P;
-        unitcell::Tuple{Int, Int} = (1, 1)
-    ) where {P <: ProductSpace}
-    return ProductStateEnv(f, T, N, fill(D_north, unitcell), fill(D_east, unitcell))
+function ProductStateEnv(Ds_north::A, args...; kwargs...) where {A <: AbstractMatrix{<:VectorSpace}}
+    return ProductStateEnv(randn, ComplexF64, Ds_north, args...; kwargs...)
 end
 
 """
-    ProductStateEnv([f=ones, T=ComplexF64], network::InfiniteSquareNetwork)
+    ProductStateEnv([f=randn, T=ComplexF64], network::InfiniteSquareNetwork)
 
 Construct a product state environment by specifying a corresponding [`InfiniteSquareNetwork`](@ref).
 """
@@ -94,35 +84,20 @@ function ProductStateEnv(f, T, network::InfiniteSquareNetwork)
     Ds_east = _east_edge_physical_spaces(network)
     return ProductStateEnv(f, T, Ds_north, Ds_east)
 end
-function ProductStateEnv(network::InfiniteSquareNetwork)
-    return ProductStateEnv(ones, scalartype(network), network) # TODO: do we want to use a different default function?
-end
-
-function ProductStateEnv(state::Union{InfinitePartitionFunction, InfinitePEPS, InfinitePEPO}, args...; kwargs...)
-    return ProductStateEnv(InfiniteSquareNetwork(state), args...; kwargs...)
-end
-function ProductStateEnv(state::Union{InfinitePEPS, InfinitePEPO}, args...; kwargs...)
-    return ProductStateEnv(InfiniteSquareNetwork(state), args...; kwargs...)
+function ProductStateEnv(network::Union{InfiniteSquareNetwork, InfinitePartitionFunction, InfinitePEPS})
+    return ProductStateEnv(randn, scalartype(network), network)
 end
-function ProductStateEnv(f, T, state::Union{InfinitePartitionFunction, InfinitePEPS, InfinitePEPO}, args...; kwargs...)
-    return ProductStateEnv(f, T, InfiniteSquareNetwork(state), args...; kwargs...)
+function ProductStateEnv(f, T, state::Union{InfinitePartitionFunction, InfinitePEPS}, args...)
+    return ProductStateEnv(f, T, InfiniteSquareNetwork(state), args...)
 end
 
 Base.eltype(::Type{ProductStateEnv{T}}) where {T} = T
 Base.size(env::ProductStateEnv, args...) = size(env.edges, args...)
 Base.getindex(env::ProductStateEnv, args...) = Base.getindex(env.edges, args...)
-Base.axes(env::ProductStateEnv, args...) = Base.axes(env.edges, args...)
-Base.eachindex(env::ProductStateEnv) = eachindex(IndexCartesian(), env.edges)
+Base.eachindex(index_style, env::ProductStateEnv) = eachindex(index_style, env.edges)
 VectorInterface.scalartype(::Type{ProductStateEnv{T}}) where {T} = scalartype(T)
 TensorKit.spacetype(::Type{ProductStateEnv{T}}) where {T} = spacetype(T)
 
-function eachcoordinate(x::ProductStateEnv)
-    return collect(Iterators.product(axes(x, 2), axes(x, 3)))
-end
-function eachcoordinate(x::ProductStateEnv, dirs)
-    return collect(Iterators.product(dirs, axes(x, 2), axes(x, 3)))
-end
-
 # conversion to CTMRGEnv
 """
     CTMRGEnv(prod_env::ProductStateEnv)
@@ -131,10 +106,10 @@ Construct a CTMRG environment with a trivial virtual space of bond dimension χ
 from the product state environment `prod_env`.
 """
 function CTMRGEnv(prod_env::ProductStateEnv)
-    edges = map(CartesianIndices(prod_env.edges)) do idx
-        return insertleftunit(insertleftunit(prod_env.edges[idx]), 1)
+    edges = map(eachindex(IndexCartesian(), prod_env)) do idx
+        return insertleftunit(insertleftunit(prod_env[idx]), 1)
     end
-    corners = map(CartesianIndices(edges)) do _
+    corners = map(eachindex(IndexCartesian(), prod_env)) do _
         return TensorKit.id(scalartype(prod_env), oneunit(spacetype(prod_env)))
     end
     return CTMRGEnv(corners, edges)
diff --git a/test/ctmrg/initialization.jl b/test/ctmrg/initialization.jl
index 0969cbf46..7fa7d1a30 100644
--- a/test/ctmrg/initialization.jl
+++ b/test/ctmrg/initialization.jl
@@ -47,6 +47,15 @@ boundary_alg = (; alg = :SimultaneousCTMRG, tol, verbosity, trunc, maxiter)
     env_appl, info = leading_boundary(env0_appl, n; boundary_alg...)
     @test info.convergence_error ≤ tol
 
+    # specific custom starting product state
+    p_data = ComplexF64[1 0]
+    p = Tensor(p_data, P)
+    prod_env0 = ProductStateEnv(reshape([p, p, flip(p, 1), flip(p, 1)], 4, 1, 1))
+    env0_custom = initialize_ctmrg_environment(n, ApplicationInitialization(), prod_env0)
+    # or just CTMRGEnv(prod_env0)
+    env_custom, info = leading_boundary(env0_custom, n; boundary_alg...)
+    @test info.convergence_error ≤ tol
+
     # PEPS-specific identity initialization; should throw when used on partition functions
     Random.seed!(sd)
     @test_throws ArgumentError env0_prod_id = initialize_ctmrg_environment(n, IdentityInitialization())

From e05f66290821ea6f895db8a0e89b281b6fce484c Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@ugent.be>
Date: Tue, 16 Jun 2026 14:35:39 +0200
Subject: [PATCH 17/24] Small cleanup for `BPEnv`

---
 src/environments/bp_environments.jl | 38 +++++++++--------------------
 1 file changed, 11 insertions(+), 27 deletions(-)

diff --git a/src/environments/bp_environments.jl b/src/environments/bp_environments.jl
index affb71e4b..bc0429e0a 100644
--- a/src/environments/bp_environments.jl
+++ b/src/environments/bp_environments.jl
@@ -51,7 +51,7 @@ end
 
 """
     BPEnv(
-        [f=randn, T=ComplexF64], Ds_north::A, Ds_east::A; posdef::Bool = true
+        [f=isomorphism, T=ComplexF64], Ds_north::A, Ds_east::A; posdef::Bool = true
     ) where {A <: AbstractMatrix{<:ProductSpace}}
 
 Construct a BP environment by specifying matrices of north and east virtual spaces of the
@@ -64,15 +64,9 @@ Each entry of the `Ds_north` and `Ds_east` matrices corresponds to an effective
 of the network, and can be represented as a `ProductSpace` (e.g.
 for the case of a network representing overlaps of PEPSs).
 """
-function BPEnv(
-        Ds_north::A, Ds_east::A; posdef::Bool = true
-    ) where {A <: AbstractMatrix{<:ProductSpace}}
-    return BPEnv(randn, ComplexF64, N, Ds_north, Ds_east; posdef)
-end
 function BPEnv(
         f, T, Ds_north::A, Ds_east::A; posdef::Bool = true
     ) where {A <: AbstractMatrix{<:ProductSpace}}
-    # no recursive broadcasting?
     Ds_south = _elementwise_dual.(circshift(Ds_north, (-1, 0)))
     Ds_west = _elementwise_dual.(circshift(Ds_east, (0, 1)))
     messages = map(Iterators.product(1:4, axes(Ds_north, 1), axes(Ds_north, 2))) do (dir, r, c)
@@ -93,10 +87,15 @@ function BPEnv(
     normalize!.(messages)
     return BPEnv(messages)
 end
+function BPEnv(
+        D_north::P, args...; kwargs...
+    ) where {P <: Union{Matrix{ProductSpace}, ProductSpace}}
+    return BPEnv(isomorphism, ComplexF64, D_north, args...; kwargs...)
+end
 
 """
     BPEnv(
-        [f=randn, T=ComplexF64], D_north::P, D_east::P;
+        [f=isomorphism, T=ComplexF64], D_north::P, D_east::P;
         unitcell::Tuple{Int, Int} = (1, 1), posdef::Bool = true
     ) where {P <: ProductSpace}
 
@@ -104,21 +103,15 @@ Construct a BP environment by specifying the north and east virtual spaces of th
 corresponding [`InfiniteSquareNetwork`](@ref). The network unit cell can be specified
 by the `unitcell` keyword argument.
 """
-function BPEnv(
-        D_north::P, D_east::P;
-        unitcell::Tuple{Int, Int} = (1, 1), posdef::Bool = true
-    ) where {P <: ProductSpace}
-    return BPEnv(randn, ComplexF64, D_north, D_east; unitcell, posdef)
-end
 function BPEnv(
         f, T, D_north::P, D_east::P;
         unitcell::Tuple{Int, Int} = (1, 1), posdef::Bool = true
     ) where {P <: ProductSpace}
-    return BPEnv(f, T, N, fill(D_north, unitcell), fill(D_east, unitcell); posdef)
+    return BPEnv(f, T, N, _fill_edge_physical_spaces(D_north, D_east; unitcell)...; posdef)
 end
 
 """
-    BPEnv([f=randn, T=ComplexF64], network::InfiniteSquareNetwork; posdef::Bool = true)
+    BPEnv([f=isomorphism, T=ComplexF64], network::InfiniteSquareNetwork; posdef::Bool = true)
 
 Construct a BP environment by specifying a corresponding [`InfiniteSquareNetwork`](@ref).
 """
@@ -127,17 +120,8 @@ function BPEnv(f, T, network::InfiniteSquareNetwork; posdef::Bool = true)
     Ds_east = _east_edge_physical_spaces(network)
     return BPEnv(f, T, Ds_north, Ds_east; posdef)
 end
-function BPEnv(network::InfiniteSquareNetwork; posdef::Bool = true)
-    return BPEnv(randn, scalartype(network), network; posdef)
-end
-
-function BPEnv(state::InfinitePartitionFunction, args...; kwargs...)
-    return BPEnv(InfiniteSquareNetwork(state), args...; kwargs...)
-end
-function BPEnv(state::Union{InfinitePEPS, InfinitePEPO}, args...; kwargs...)
-    bp_env = BPEnv(InfiniteSquareNetwork(state), args...; kwargs...)
-    TensorKit.id!.(bp_env.messages)
-    return bp_env
+function BPEnv(network::Union{InfiniteSquareNetwork, InfinitePartitionFunction, InfinitePEPS, InfinitePEPO}, args...; kwargs...)
+    return BPEnv(isomorphism, scalartype(network), network, args...; kwargs...)
 end
 function BPEnv(f, T, state::Union{InfinitePartitionFunction, InfinitePEPS, InfinitePEPO}, args...; kwargs...)
     return BPEnv(f, T, InfiniteSquareNetwork(state), args...; kwargs...)

From af90a1345b02bc81fc53cc28e998d0acb1ad08c4 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@ugent.be>
Date: Wed, 17 Jun 2026 10:56:08 +0200
Subject: [PATCH 18/24] TensorMap data should be a matrix

---
 test/ctmrg/initialization.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/ctmrg/initialization.jl b/test/ctmrg/initialization.jl
index 7fa7d1a30..98d0d6e7b 100644
--- a/test/ctmrg/initialization.jl
+++ b/test/ctmrg/initialization.jl
@@ -48,7 +48,7 @@ boundary_alg = (; alg = :SimultaneousCTMRG, tol, verbosity, trunc, maxiter)
     @test info.convergence_error ≤ tol
 
     # specific custom starting product state
-    p_data = ComplexF64[1 0]
+    p_data = ComplexF64[1; 0;;]
     p = Tensor(p_data, P)
     prod_env0 = ProductStateEnv(reshape([p, p, flip(p, 1), flip(p, 1)], 4, 1, 1))
     env0_custom = initialize_ctmrg_environment(n, ApplicationInitialization(), prod_env0)

From 9dfd6b302735ccf492dc7572568502e7d3c3645f Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@ugent.be>
Date: Wed, 17 Jun 2026 11:53:09 +0200
Subject: [PATCH 19/24] Forgot import

---
 test/ctmrg/initialization.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/ctmrg/initialization.jl b/test/ctmrg/initialization.jl
index 98d0d6e7b..78d6211dc 100644
--- a/test/ctmrg/initialization.jl
+++ b/test/ctmrg/initialization.jl
@@ -4,6 +4,7 @@ using PEPSKit
 using Random
 
 using MPSKitModels: classical_ising
+using PEPSKit: ProductStateEnv
 
 sd = 12345
 

From d0ec8d95e6735a686fa70392d2507fb7eb8b2913 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@ugent.be>
Date: Wed, 17 Jun 2026 13:33:51 +0200
Subject: [PATCH 20/24] Actually use seed

---
 test/ctmrg/initialization.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/ctmrg/initialization.jl b/test/ctmrg/initialization.jl
index 78d6211dc..577fb30f5 100644
--- a/test/ctmrg/initialization.jl
+++ b/test/ctmrg/initialization.jl
@@ -64,6 +64,7 @@ end
 
 @testset "CTMRG environment initialization for PEPS with $S symmetry" for S in symmetries
     # initialize
+    Random.seed!(sd)
     P = make_space(S, d)
     Vpeps = make_space(S, D)
     Venv = make_space(S, χ)

From 9a72484b00f2a05d6cc3cb2303cbd24283093ced Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Fri, 19 Jun 2026 12:56:10 +0200
Subject: [PATCH 21/24] Only seed at the start

---
 test/ctmrg/initialization.jl | 9 +--------
 1 file changed, 1 insertion(+), 8 deletions(-)

diff --git a/test/ctmrg/initialization.jl b/test/ctmrg/initialization.jl
index 577fb30f5..fcb4ae7b0 100644
--- a/test/ctmrg/initialization.jl
+++ b/test/ctmrg/initialization.jl
@@ -24,6 +24,7 @@ boundary_alg = (; alg = :SimultaneousCTMRG, tol, verbosity, trunc, maxiter)
 
 @testset "CTMRG environment initialization for critical ising with $S symmetry (#255)" for S in symmetries
     # initialize
+    Random.seed!(sd)
     T = classical_ising(S)
     O = T[1]
     n = InfinitePartitionFunction([O O; O O])
@@ -31,19 +32,16 @@ boundary_alg = (; alg = :SimultaneousCTMRG, tol, verbosity, trunc, maxiter)
     P = space(O, 2)
 
     # random, doesn't converge
-    Random.seed!(sd)
     env0_rand = initialize_ctmrg_environment(n, RandomInitialization(), Venv)
     env_rand, info = leading_boundary(env0_rand, n; boundary_alg...)
     @test_broken info.convergence_error ≤ tol
 
     # embedded random product state, converges
-    Random.seed!(sd)
     env0_prod = initialize_ctmrg_environment(n, ProductStateInitialization())
     env_prod, info = leading_boundary(env0_prod, n; boundary_alg...)
     @test info.convergence_error ≤ tol
 
     # grown product state, converges
-    Random.seed!(sd)
     env0_appl = initialize_ctmrg_environment(n, ApplicationInitialization())
     env_appl, info = leading_boundary(env0_appl, n; boundary_alg...)
     @test info.convergence_error ≤ tol
@@ -58,7 +56,6 @@ boundary_alg = (; alg = :SimultaneousCTMRG, tol, verbosity, trunc, maxiter)
     @test info.convergence_error ≤ tol
 
     # PEPS-specific identity initialization; should throw when used on partition functions
-    Random.seed!(sd)
     @test_throws ArgumentError env0_prod_id = initialize_ctmrg_environment(n, IdentityInitialization())
 end
 
@@ -72,25 +69,21 @@ end
     n = InfiniteSquareNetwork(peps)
 
     # random, converges
-    Random.seed!(sd)
     env0_rand = initialize_ctmrg_environment(n, RandomInitialization(), Venv)
     env_rand, info = leading_boundary(env0_rand, n; boundary_alg...)
     @test info.convergence_error ≤ tol
 
     # embedded random product state, converges
-    Random.seed!(sd)
     env0_prod = initialize_ctmrg_environment(n, ProductStateInitialization())
     env_prod, info = leading_boundary(env0_prod, n; boundary_alg...)
     @test info.convergence_error ≤ tol
 
     # embedded product state as identity from ket to bra, converges
-    Random.seed!(sd)
     env0_prod_id = initialize_ctmrg_environment(n, IdentityInitialization())
     env_prod, info = leading_boundary(env0_prod_id, n; boundary_alg...)
     @test info.convergence_error ≤ tol
 
     # grown product state, converges
-    Random.seed!(sd)
     env0_appl = initialize_ctmrg_environment(n, ApplicationInitialization())
     env_appl, info = leading_boundary(env0_appl, n; boundary_alg...)
     @test info.convergence_error ≤ tol

From f51ed73fbd32ea73f76c3a892a5d7c8d15ea4ac5 Mon Sep 17 00:00:00 2001
From: Lander Burgelman <39218680+leburgel@users.noreply.github.com>
Date: Tue, 23 Jun 2026 13:30:42 +0200
Subject: [PATCH 22/24] Update src/environments/product_state_environments.jl

Co-authored-by: Lukas Devos <ldevos98@gmail.com>
---
 src/environments/product_state_environments.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/environments/product_state_environments.jl b/src/environments/product_state_environments.jl
index fd660d659..ea6550e10 100644
--- a/src/environments/product_state_environments.jl
+++ b/src/environments/product_state_environments.jl
@@ -110,7 +110,7 @@ function CTMRGEnv(prod_env::ProductStateEnv)
         return insertleftunit(insertleftunit(prod_env[idx]), 1)
     end
     corners = map(eachindex(IndexCartesian(), prod_env)) do _
-        return TensorKit.id(scalartype(prod_env), oneunit(spacetype(prod_env)))
+        return TensorKit.id(storagetype(prod_env), oneunit(spacetype(prod_env)))
     end
     return CTMRGEnv(corners, edges)
 end

From b8a36ddbe64dbbcc2eea33fba7dbb300cca496b0 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Tue, 23 Jun 2026 17:22:42 +0200
Subject: [PATCH 23/24] Define `TensorKit.storagetype` for `ProductStateEnv`

---
 src/environments/product_state_environments.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/environments/product_state_environments.jl b/src/environments/product_state_environments.jl
index ea6550e10..394bc5751 100644
--- a/src/environments/product_state_environments.jl
+++ b/src/environments/product_state_environments.jl
@@ -96,6 +96,7 @@ Base.size(env::ProductStateEnv, args...) = size(env.edges, args...)
 Base.getindex(env::ProductStateEnv, args...) = Base.getindex(env.edges, args...)
 Base.eachindex(index_style, env::ProductStateEnv) = eachindex(index_style, env.edges)
 VectorInterface.scalartype(::Type{ProductStateEnv{T}}) where {T} = scalartype(T)
+TensorKit.storagetype(::Type{ProductStateEnv{T}}) where {T} = storagetype(T)
 TensorKit.spacetype(::Type{ProductStateEnv{T}}) where {T} = spacetype(T)
 
 # conversion to CTMRGEnv

From bed9c2f80dae3cdeeaa0154179e825517df6f8db Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Wed, 24 Jun 2026 09:45:51 +0200
Subject: [PATCH 24/24] Add actual `CTMRGEnv(::CTMRGEnv)` constructor

---
 src/algorithms/ctmrg/initialization.jl | 5 +----
 src/environments/ctmrg_environments.jl | 3 +++
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/algorithms/ctmrg/initialization.jl b/src/algorithms/ctmrg/initialization.jl
index 73e7b0769..235611e0b 100644
--- a/src/algorithms/ctmrg/initialization.jl
+++ b/src/algorithms/ctmrg/initialization.jl
@@ -28,9 +28,6 @@ function initialize_ctmrg_environment(
     return env
 end
 
-_CTMRGEnv(env) = CTMRGEnv(env)
-_CTMRGEnv(env::CTMRGEnv) = env
-
 """
     initialize_ctmrg_environment([elt::Type{<:Number},] n::InfiniteSquareNetwork, alg::RandomInitialization, [env0])
 
@@ -45,7 +42,7 @@ function initialize_ctmrg_environment(
         env0 = ProductStateEnv(alg.f, elt, n)
     )
     dummy_alg = SimultaneousCTMRG(trunc = (; alg = :notrunc))
-    env, = ctmrg_iteration(n, _CTMRGEnv(env0), dummy_alg)
+    env, = ctmrg_iteration(n, CTMRGEnv(env0), dummy_alg)
     return env
 end
 
diff --git a/src/environments/ctmrg_environments.jl b/src/environments/ctmrg_environments.jl
index 17b4594d1..cabc0e764 100644
--- a/src/environments/ctmrg_environments.jl
+++ b/src/environments/ctmrg_environments.jl
@@ -230,6 +230,9 @@ function CTMRGEnv(f, T, state::Union{InfinitePartitionFunction, InfinitePEPS}, a
     return CTMRGEnv(f, T, InfiniteSquareNetwork(state), args...)
 end
 
+# copy-like constructor
+CTMRGEnv(env::CTMRGEnv) = CTMRGEnv(env.corners, env.edges)
+
 @non_differentiable CTMRGEnv(state::Union{InfinitePartitionFunction, InfinitePEPS}, args...)
 
 # Custom adjoint for CTMRGEnv constructor, needed for fixed-point differentiation