diff --git a/README.md b/README.md
index 381cb8d90..fee1a4f5a 100644
--- a/README.md
+++ b/README.md
@@ -54,9 +54,9 @@ opt_alg = PEPSOptimize(;
 
 # ground state search
 state = InfinitePEPS(2, D)
-ctm = leading_boundary(CTMRGEnv(state, ComplexSpace(chi)), state, ctm_alg)
-result = fixedpoint(state, H, opt_alg, ctm)
+ctm, = leading_boundary(CTMRGEnv(state, ComplexSpace(chi)), state, ctm_alg)
+peps, env, E, = fixedpoint(H, state, ctm, opt_alg)
 
-@show result.E # -0.6625...
+@show E # -0.6625...
 ```
 
diff --git a/docs/src/index.md b/docs/src/index.md
index 115349e12..dd1d66e59 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -36,8 +36,8 @@ opt_alg = PEPSOptimize(;
 
 # ground state search
 state = InfinitePEPS(2, D)
-ctm = leading_boundary(CTMRGEnv(state, ComplexSpace(chi)), state, ctm_alg)
-result = fixedpoint(state, H, opt_alg, ctm)
+ctm, = leading_boundary(CTMRGEnv(state, ComplexSpace(chi)), state, ctm_alg)
+peps, env, E, = fixedpoint(H, state, ctm, opt_alg)
 
-@show result.E # -0.6625...
+@show E # -0.6625...
 ```
diff --git a/examples/boundary_mps.jl b/examples/boundary_mps.jl
index 3047d606e..1f8ccbbb3 100644
--- a/examples/boundary_mps.jl
+++ b/examples/boundary_mps.jl
@@ -26,13 +26,13 @@ T = InfiniteTransferPEPS(peps, 1, 1)
 mps = PEPSKit.initializeMPS(T, [ComplexSpace(20)])
 
 # We then find the leading boundary MPS fixed point using the VUMPS algorithm
-mps, envs, ϵ = leading_boundary(mps, T, VUMPS())
+mps, env, ϵ = leading_boundary(mps, T, VUMPS())
 
 # The norm of the state per unit cell is then given by the expectation value <mps|T|mps>
 N = abs(prod(expectation_value(mps, T)))
 
 # This can be compared to the result obtained using the CTMRG algorithm
-ctm = leading_boundary(
+ctm, = leading_boundary(
     peps, SimultaneousCTMRG(; verbosity=1), CTMRGEnv(peps, ComplexSpace(20))
 )
 N´ = abs(norm(peps, ctm))
@@ -52,10 +52,10 @@ peps2 = InfinitePEPS(ComplexSpace(2), ComplexSpace(2); unitcell=(2, 2))
 T2 = PEPSKit.MultilineTransferPEPS(peps2, 1)
 
 mps2 = PEPSKit.initializeMPS(T2, fill(ComplexSpace(20), 2, 2))
-mps2, envs2, ϵ = leading_boundary(mps2, T2, VUMPS())
+mps2, env2, ϵ = leading_boundary(mps2, T2, VUMPS())
 N2 = abs(prod(expectation_value(mps2, T2)))
 
-ctm2 = leading_boundary(
+ctm2, = leading_boundary(
     peps2, SimultaneousCTMRG(; verbosity=1), CTMRGEnv(peps2, ComplexSpace(20))
 )
 N2´ = abs(norm(peps2, ctm2))
@@ -77,7 +77,7 @@ pepo = ising_pepo(1)
 T3 = InfiniteTransferPEPO(peps, pepo, 1, 1)
 
 mps3 = PEPSKit.initializeMPS(T3, [ComplexSpace(20)])
-mps3, envs3, ϵ = leading_boundary(mps3, T3, VUMPS())
+mps3, env3, ϵ = leading_boundary(mps3, T3, VUMPS())
 @show N3 = abs(prod(expectation_value(mps3, T3)))
 
 # These objects and routines can be used to optimize PEPS fixed points of 3D partition
diff --git a/examples/heisenberg.jl b/examples/heisenberg.jl
index 3b0b0d657..d54a7eb39 100644
--- a/examples/heisenberg.jl
+++ b/examples/heisenberg.jl
@@ -26,6 +26,6 @@ opt_alg = PEPSOptimize(;
 # E/N = −0.6694421, which is a QMC estimate from https://arxiv.org/abs/1101.3281.
 # Of course there is a noticable bias for small χbond and χenv.
 ψ₀ = InfinitePEPS(2, χbond)
-env₀ = leading_boundary(CTMRGEnv(ψ₀, ℂ^χenv), ψ₀, ctm_alg)
-result = fixedpoint(ψ₀, H, opt_alg, env₀)
-@show result.E
+env₀, = leading_boundary(CTMRGEnv(ψ₀, ℂ^χenv), ψ₀, ctm_alg)
+peps, env, E, = fixedpoint(H, ψ, env₀, opt_alg₀)
+@show E
diff --git a/examples/heisenberg_evol/heis_tools.jl b/examples/heisenberg_evol/heis_tools.jl
index 80aec6f9a..df68cee4b 100644
--- a/examples/heisenberg_evol/heis_tools.jl
+++ b/examples/heisenberg_evol/heis_tools.jl
@@ -16,7 +16,7 @@ using PEPSKit
 """
 Measure magnetization on each site
 """
-function cal_mags(peps::InfinitePEPS, envs::CTMRGEnv)
+function cal_mags(peps::InfinitePEPS, env::CTMRGEnv)
     N1, N2 = size(peps)
     lattice = collect(space(t, 1) for t in peps.A)
     # detect symmetry on physical axis
@@ -32,7 +32,7 @@ function cal_mags(peps::InfinitePEPS, envs::CTMRGEnv)
     return [
         collect(
             expectation_value(
-                peps, LocalOperator(lattice, (CartesianIndex(r, c),) => Sa), envs
+                peps, LocalOperator(lattice, (CartesianIndex(r, c),) => Sa), env
             ) for (r, c) in Iterators.product(1:N1, 1:N2)
         ) for Sa in Sas
     ]
@@ -41,11 +41,11 @@ end
 """
 Measure physical quantities for Heisenberg model
 """
-function measure_heis(peps::InfinitePEPS, H::LocalOperator, envs::CTMRGEnv)
+function measure_heis(peps::InfinitePEPS, H::LocalOperator, env::CTMRGEnv)
     results = Dict{String,Any}()
     N1, N2 = size(peps)
-    results["e_site"] = costfun(peps, envs, H) / (N1 * N2)
-    results["mag"] = cal_mags(peps, envs)
+    results["e_site"] = costfun(peps, env, H) / (N1 * N2)
+    results["mag"] = cal_mags(peps, env)
     results["mag_norm"] = collect(
         norm([mags[r, c] for mags in results["mag"]]) for
         (r, c) in Iterators.product(1:N1, 1:N2)
diff --git a/examples/heisenberg_evol/simpleupdate.jl b/examples/heisenberg_evol/simpleupdate.jl
index b16e6868b..ee8fedba2 100644
--- a/examples/heisenberg_evol/simpleupdate.jl
+++ b/examples/heisenberg_evol/simpleupdate.jl
@@ -3,7 +3,7 @@ include("heis_tools.jl")
 # benchmark data is from Phys. Rev. B 94, 035133 (2016)
 Dbond, χenv, symm = 4, 16, Trivial
 trscheme_peps = truncerr(1e-10) & truncdim(Dbond)
-trscheme_envs = truncerr(1e-10) & truncdim(χenv)
+trscheme_env = truncerr(1e-10) & truncdim(χenv)
 N1, N2 = 2, 2
 # Heisenberg model Hamiltonian
 # (already only includes nearest neighbor terms)
@@ -41,10 +41,10 @@ for (dt, tol) in zip(dts, tols)
 end
 # measure physical quantities with CTMRG
 peps_ = InfinitePEPS(peps)
-envs = CTMRGEnv(rand, Float64, peps_, Espace)
-ctm_alg = SequentialCTMRG(; tol=1e-10, verbosity=2, trscheme=trscheme_envs)
-envs = leading_boundary(envs, peps_, ctm_alg)
-meas = measure_heis(peps_, ham, envs)
+env₀ = CTMRGEnv(rand, Float64, peps_, Espace)
+ctm_alg = SequentialCTMRG(; tol=1e-10, verbosity=2, trscheme=trscheme_env)
+env = leading_boundary(env₀, peps_, ctm_alg)
+meas = measure_heis(peps_, ham, env)
 display(meas)
 @info @sprintf("Energy = %.8f\n", meas["e_site"])
 @info @sprintf("Staggered magnetization = %.8f\n", mean(meas["mag_norm"]))
diff --git a/examples/hubbard_su.jl b/examples/hubbard_su.jl
index 6478a28e0..9c1468be5 100644
--- a/examples/hubbard_su.jl
+++ b/examples/hubbard_su.jl
@@ -40,10 +40,10 @@ peps = InfinitePEPS(peps)
 # CTMRG
 χenv0, χenv = 6, 20
 Espace = Vect[fℤ₂](0 => χenv0 / 2, 1 => χenv0 / 2)
-envs = CTMRGEnv(randn, Float64, peps, Espace)
+env = CTMRGEnv(randn, Float64, peps, Espace)
 for χ in [χenv0, χenv]
     ctm_alg = SequentialCTMRG(; maxiter=300, tol=1e-7)
-    envs = leading_boundary(envs, peps, ctm_alg)
+    env, = leading_boundary(env, peps, ctm_alg)
 end
 
 # Benchmark values of the ground state energy from
@@ -53,7 +53,7 @@ Es_exact = Dict(0 => -1.62, 2 => -0.176, 4 => 0.8603, 6 => -0.6567, 8 => -0.5243
 E_exact = Es_exact[U] - U / 2
 
 # measure energy
-E = costfun(peps, envs, ham) / (N1 * N2)
+E = cost_function(peps, env, ham) / (N1 * N2)
 @info "Energy           = $E"
 @info "Benchmark energy = $E_exact"
 @test isapprox(E, E_exact; atol=5e-2)
diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 06dbf76f0..caae82240 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -54,56 +54,65 @@ include("algorithms/time_evolution/simpleupdate.jl")
 
 include("algorithms/toolbox.jl")
 
-include("algorithms/peps_opt.jl")
-
 include("utility/symmetrization.jl")
 
+include("algorithms/optimization/fixed_point_differentiation.jl")
+include("algorithms/optimization/peps_optimization.jl")
+
 """
     module Defaults
-        const ctmrg_maxiter = 100
-        const ctmrg_miniter = 4
-        const ctmrg_tol = 1e-8
-        const fpgrad_maxiter = 30
-        const fpgrad_tol = 1e-6
-        const reuse_env = true
-        const trscheme = FixedSpaceTruncation()
-        const fwd_alg = TensorKit.SDD()
-        const rrule_alg = Arnoldi(; tol=1e-2fpgrad_tol, krylovdim=48, verbosity=-1)
-        const svd_alg = SVDAdjoint(; fwd_alg, rrule_alg)
-        const projector_alg_type = HalfInfiniteProjector
-        const projector_alg = projector_alg_type(svd_alg, trscheme, 2)
-        const ctmrg_alg = SimultaneousCTMRG(
-            ctmrg_tol, ctmrg_maxiter, ctmrg_miniter, 2, projector_alg
-        )
-        const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=3)
-        const gradient_linsolver = KrylovKit.BiCGStab(;
-            maxiter=Defaults.fpgrad_maxiter, tol=Defaults.fpgrad_tol
-        )
-        const iterscheme = :fixed
-        const gradient_alg = LinSolver(; solver=gradient_linsolver, iterscheme)
-        const scheduler = Ref{Scheduler}(Threads.nthreads() == 1 ? SerialScheduler() : DynamicScheduler())
-    end
 
-Module containing default values that represent typical algorithm parameters.
+Module containing default algorithm parameter values and arguments.
 
-- `ctmrg_maxiter`: Maximal number of CTMRG iterations per run
-- `ctmrg_miniter`: Minimal number of CTMRG carried out
-- `ctmrg_tol`: Tolerance checking singular value and norm convergence
-- `fpgrad_maxiter`: Maximal number of iterations for computing the CTMRG fixed-point gradient
-- `fpgrad_tol`: Convergence tolerance for the fixed-point gradient iteration
-- `reuse_env`: If `true`, the current optimization step is initialized on the previous environment
-- `trscheme`: Truncation scheme for SVDs and other decompositions
-- `fwd_alg`: SVD algorithm that is used in the forward pass
+# CTMRG
+- `ctmrg_tol=1e-8`: Tolerance checking singular value and norm convergence
+- `ctmrg_maxiter=100`: Maximal number of CTMRG iterations per run
+- `ctmrg_miniter=4`: Minimal number of CTMRG carried out
+- `trscheme=FixedSpaceTruncation()`: Truncation scheme for SVDs and other decompositions
+- `fwd_alg=TensorKit.SDD()`: SVD algorithm that is used in the forward pass
 - `rrule_alg`: Reverse-rule for differentiating that SVD
-- `svd_alg`: Combination of `fwd_alg` and `rrule_alg`
-- `projector_alg_type`: Default type of projector algorithm
+
+    ```
+    rrule_alg = Arnoldi(; tol=ctmrg_tol, krylovdim=48, verbosity=-1)
+    ```
+
+- `svd_alg=SVDAdjoint(; fwd_alg, rrule_alg)`: Combination of `fwd_alg` and `rrule_alg`
+- `projector_alg_type=HalfInfiniteProjector`: Default type of projector algorithm
 - `projector_alg`: Algorithm to compute CTMRG projectors
+
+    ```
+    projector_alg = projector_alg_type(; svd_alg, trscheme, verbosity=0)
+    ```
+
 - `ctmrg_alg`: Algorithm for performing CTMRG runs
-- `optimizer`: Optimization algorithm for PEPS ground-state optimization
+
+    ```
+    ctmrg_alg = SimultaneousCTMRG(
+        ctmrg_tol, ctmrg_maxiter, ctmrg_miniter, 2, projector_alg
+    )
+    ```
+
+# Optimization
+- `fpgrad_maxiter=30`: Maximal number of iterations for computing the CTMRG fixed-point gradient
+- `fpgrad_tol=1e-6`: Convergence tolerance for the fixed-point gradient iteration
+- `iterscheme=:fixed`: Scheme for differentiating one CTMRG iteration
 - `gradient_linsolver`: Default linear solver for the `LinSolver` gradient algorithm
-- `iterscheme`: Scheme for differentiating one CTMRG iteration
+
+    ```
+    gradient_linsolver=KrylovKit.BiCGStab(; maxiter=fpgrad_maxiter, tol=fpgrad_tol)
+    ```
+
 - `gradient_alg`: Algorithm to compute the gradient fixed-point
-- `scheduler`: Multi-threading scheduler which can be accessed via `set_scheduler!`
+
+    ```
+    gradient_alg = LinSolver(; solver=gradient_linsolver, iterscheme)
+    ```
+
+- `reuse_env=true`: If `true`, the current optimization step is initialized on the previous environment
+- `optimizer=LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=3)`: Default `OptimKit.OptimizerAlgorithm` for PEPS optimization
+
+# OhMyThreads scheduler
+- `scheduler=Ref{Scheduler}(...)`: Multi-threading scheduler which can be accessed via `set_scheduler!`
 """
 module Defaults
     using TensorKit, KrylovKit, OptimKit, OhMyThreads
@@ -113,28 +122,30 @@ module Defaults
         SVDAdjoint,
         HalfInfiniteProjector,
         SimultaneousCTMRG
+
+    # CTMRG
     const ctmrg_tol = 1e-8
     const ctmrg_maxiter = 100
     const ctmrg_miniter = 4
-    const fpgrad_maxiter = 30
-    const fpgrad_tol = 1e-6
     const sparse = false
-    const reuse_env = true
     const trscheme = FixedSpaceTruncation()
     const fwd_alg = TensorKit.SDD()
-    const rrule_alg = Arnoldi(; tol=1e-2fpgrad_tol, krylovdim=48, verbosity=-1)
+    const rrule_alg = Arnoldi(; tol=ctmrg_tol, krylovdim=48, verbosity=-1)
     const svd_alg = SVDAdjoint(; fwd_alg, rrule_alg)
     const projector_alg_type = HalfInfiniteProjector
-    const projector_alg = projector_alg_type(svd_alg, trscheme, 2)
+    const projector_alg = projector_alg_type(; svd_alg, trscheme, verbosity=0)
     const ctmrg_alg = SimultaneousCTMRG(
         ctmrg_tol, ctmrg_maxiter, ctmrg_miniter, 2, projector_alg
     )
-    const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=3)
-    const gradient_linsolver = KrylovKit.BiCGStab(;
-        maxiter=Defaults.fpgrad_maxiter, tol=Defaults.fpgrad_tol
-    )
+
+    # Optimization
+    const fpgrad_maxiter = 30
+    const fpgrad_tol = 1e-6
+    const gradient_linsolver = KrylovKit.BiCGStab(; maxiter=fpgrad_maxiter, tol=fpgrad_tol)
     const iterscheme = :fixed
     const gradient_alg = LinSolver(; solver=gradient_linsolver, iterscheme)
+    const reuse_env = true
+    const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=3)
 
     # OhMyThreads scheduler defaults
     const scheduler = Ref{Scheduler}()
@@ -187,7 +198,7 @@ export SVDAdjoint, IterSVD
 export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG
 export FixedSpaceTruncation, HalfInfiniteProjector, FullInfiniteProjector
 export LocalOperator
-export expectation_value, costfun, product_peps, correlation_length
+export expectation_value, cost_function, product_peps, correlation_length
 export leading_boundary
 export PEPSOptimize, GeomSum, ManualIter, LinSolver
 export fixedpoint
diff --git a/src/algorithms/contractions/ctmrg_contractions.jl b/src/algorithms/contractions/ctmrg_contractions.jl
index 6851babfb..2ec63fcfe 100644
--- a/src/algorithms/contractions/ctmrg_contractions.jl
+++ b/src/algorithms/contractions/ctmrg_contractions.jl
@@ -7,8 +7,8 @@ const CTMRGCornerTensor{T,S} = AbstractTensorMap{T,S,1,1}
 # ----------------------------
 
 """
-    enlarge_northwest_corner((row, col), envs, ket::InfinitePEPS, bra::InfinitePEPS=ket)
-    enlarge_northwest_corner((row, col), envs, partfunc::InfinitePartitionFunction)
+    enlarge_northwest_corner((row, col), env, ket::InfinitePEPS, bra::InfinitePEPS=ket)
+    enlarge_northwest_corner((row, col), env, partfunc::InfinitePartitionFunction)
     enlarge_northwest_corner(E_west, C_northwest, E_north, ket::PEPSTensor, bra::PEPSTensor=ket)
     enlarge_northwest_corner(E_west, C_northwest, E_north, partfunc::PartitionFunctionTensor)
 
@@ -28,11 +28,11 @@ The networks and tensors (denoted `A` below) can correspond to either:
 ```
 """
 function enlarge_northwest_corner(
-    (row, col), envs::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket
+    (row, col), env::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket
 )
-    E_west = envs.edges[WEST, row, _prev(col, end)]
-    C_northwest = envs.corners[NORTHWEST, _prev(row, end), _prev(col, end)]
-    E_north = envs.edges[NORTH, _prev(row, end), col]
+    E_west = env.edges[WEST, row, _prev(col, end)]
+    C_northwest = env.corners[NORTHWEST, _prev(row, end), _prev(col, end)]
+    E_north = env.edges[NORTH, _prev(row, end), col]
     return enlarge_northwest_corner(
         E_west, C_northwest, E_north, ket[row, col], bra[row, col]
     )
@@ -51,10 +51,10 @@ function enlarge_northwest_corner(
         ket[d; D3 D_Eabove D_Sabove D1] *
         conj(bra[d; D4 D_Ebelow D_Sbelow D2])
 end
-function enlarge_northwest_corner((row, col), envs::CTMRGEnv, partfunc::InfinitePF)
-    E_west = envs.edges[WEST, row, _prev(col, end)]
-    C_northwest = envs.corners[NORTHWEST, _prev(row, end), _prev(col, end)]
-    E_north = envs.edges[NORTH, _prev(row, end), col]
+function enlarge_northwest_corner((row, col), env::CTMRGEnv, partfunc::InfinitePF)
+    E_west = env.edges[WEST, row, _prev(col, end)]
+    C_northwest = env.corners[NORTHWEST, _prev(row, end), _prev(col, end)]
+    E_north = env.edges[NORTH, _prev(row, end), col]
     return enlarge_northwest_corner(E_west, C_northwest, E_north, partfunc[row, col])
 end
 function enlarge_northwest_corner(
@@ -71,8 +71,8 @@ function enlarge_northwest_corner(
 end
 
 """
-    enlarge_northeast_corner((row, col), envs, ket::InfinitePEPS, bra::InfinitePEPS=ket)
-    enlarge_northeast_corner((row, col), envs, partfunc::InfinitePartitionFunction)
+    enlarge_northeast_corner((row, col), env, ket::InfinitePEPS, bra::InfinitePEPS=ket)
+    enlarge_northeast_corner((row, col), env, partfunc::InfinitePartitionFunction)
     enlarge_northeast_corner(E_north, C_northeast, E_east, ket::PEPSTensor, bra::PEPSTensor=ket)
     enlarge_northeast_corner(E_north, C_northeast, E_east, partfunc::PartitionFunctionTensor)
 
@@ -92,11 +92,11 @@ The networks and tensors (denoted `A` below) can correspond to either:
 ```
 """
 function enlarge_northeast_corner(
-    (row, col), envs::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket
+    (row, col), env::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket
 )
-    E_north = envs.edges[NORTH, _prev(row, end), col]
-    C_northeast = envs.corners[NORTHEAST, _prev(row, end), _next(col, end)]
-    E_east = envs.edges[EAST, row, _next(col, end)]
+    E_north = env.edges[NORTH, _prev(row, end), col]
+    C_northeast = env.corners[NORTHEAST, _prev(row, end), _next(col, end)]
+    E_east = env.edges[EAST, row, _next(col, end)]
     return enlarge_northeast_corner(
         E_north, C_northeast, E_east, ket[row, col], bra[row, col]
     )
@@ -115,10 +115,10 @@ function enlarge_northeast_corner(
         ket[d; D1 D3 D_Sabove D_Wabove] *
         conj(bra[d; D2 D4 D_Sbelow D_Wbelow])
 end
-function enlarge_northeast_corner((row, col), envs::CTMRGEnv, partfunc::InfinitePF)
-    E_north = envs.edges[NORTH, _prev(row, end), col]
-    C_northeast = envs.corners[NORTHEAST, _prev(row, end), _next(col, end)]
-    E_east = envs.edges[EAST, row, _next(col, end)]
+function enlarge_northeast_corner((row, col), env::CTMRGEnv, partfunc::InfinitePF)
+    E_north = env.edges[NORTH, _prev(row, end), col]
+    C_northeast = env.corners[NORTHEAST, _prev(row, end), _next(col, end)]
+    E_east = env.edges[EAST, row, _next(col, end)]
     return enlarge_northeast_corner(E_north, C_northeast, E_east, partfunc[row, col])
 end
 function enlarge_northeast_corner(
@@ -135,8 +135,8 @@ function enlarge_northeast_corner(
 end
 
 """
-    enlarge_southeast_corner((row, col), envs, ket::InfinitePEPS, bra::InfinitePEPS=ket)
-    enlarge_southeast_corner((row, col), envs, partfunc::InfinitePartitionFunction)
+    enlarge_southeast_corner((row, col), env, ket::InfinitePEPS, bra::InfinitePEPS=ket)
+    enlarge_southeast_corner((row, col), env, partfunc::InfinitePartitionFunction)
     enlarge_southeast_corner(E_east, C_southeast, E_south, ket::PEPSTensor, bra::PEPSTensor=ket)
     enlarge_southeast_corner(E_east, C_southeast, E_south, partfunc::PartitionFunctionTensor)
 
@@ -156,11 +156,11 @@ The networks and tensors (denoted `A` below) can correspond to either:
 ```
 """
 function enlarge_southeast_corner(
-    (row, col), envs::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket
+    (row, col), env::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket
 )
-    E_east = envs.edges[EAST, row, _next(col, end)]
-    C_southeast = envs.corners[SOUTHEAST, _next(row, end), _next(col, end)]
-    E_south = envs.edges[SOUTH, _next(row, end), col]
+    E_east = env.edges[EAST, row, _next(col, end)]
+    C_southeast = env.corners[SOUTHEAST, _next(row, end), _next(col, end)]
+    E_south = env.edges[SOUTH, _next(row, end), col]
     return enlarge_southeast_corner(
         E_east, C_southeast, E_south, ket[row, col], bra[row, col]
     )
@@ -179,10 +179,10 @@ function enlarge_southeast_corner(
         ket[d; D_Nabove D1 D3 D_Wabove] *
         conj(bra[d; D_Nbelow D2 D4 D_Wbelow])
 end
-function enlarge_southeast_corner((row, col), envs::CTMRGEnv, partfunc::InfinitePF)
-    E_east = envs.edges[EAST, row, _next(col, end)]
-    C_southeast = envs.corners[SOUTHEAST, _next(row, end), _next(col, end)]
-    E_south = envs.edges[SOUTH, _next(row, end), col]
+function enlarge_southeast_corner((row, col), env::CTMRGEnv, partfunc::InfinitePF)
+    E_east = env.edges[EAST, row, _next(col, end)]
+    C_southeast = env.corners[SOUTHEAST, _next(row, end), _next(col, end)]
+    E_south = env.edges[SOUTH, _next(row, end), col]
     return enlarge_southeast_corner(E_east, C_southeast, E_south, partfunc[row, col])
 end
 function enlarge_southeast_corner(
@@ -199,8 +199,8 @@ function enlarge_southeast_corner(
 end
 
 """
-    enlarge_southwest_corner((row, col), envs, ket::InfinitePEPS, bra::InfinitePEPS=ket)
-    enlarge_southwest_corner((row, col), envs, partfunc::InfinitePartitionFunction)
+    enlarge_southwest_corner((row, col), env, ket::InfinitePEPS, bra::InfinitePEPS=ket)
+    enlarge_southwest_corner((row, col), env, partfunc::InfinitePartitionFunction)
     enlarge_southwest_corner(E_south, C_southwest, E_west, ket::PEPSTensor, bra::PEPSTensor=ket)
     enlarge_southwest_corner(E_south, C_southwest, E_west, partfunc::PartitionFunctionTensor)
 
@@ -220,11 +220,11 @@ The networks and tensors (denoted `A` below) can correspond to either:
 ```
 """
 function enlarge_southwest_corner(
-    (row, col), envs::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket
+    (row, col), env::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket
 )
-    E_south = envs.edges[SOUTH, _next(row, end), col]
-    C_southwest = envs.corners[SOUTHWEST, _next(row, end), _prev(col, end)]
-    E_west = envs.edges[WEST, row, _prev(col, end)]
+    E_south = env.edges[SOUTH, _next(row, end), col]
+    C_southwest = env.corners[SOUTHWEST, _next(row, end), _prev(col, end)]
+    E_west = env.edges[WEST, row, _prev(col, end)]
     return enlarge_southwest_corner(
         E_south, C_southwest, E_west, ket[row, col], bra[row, col]
     )
@@ -243,10 +243,10 @@ function enlarge_southwest_corner(
         ket[d; D_Nabove D_Eabove D1 D3] *
         conj(bra[d; D_Nbelow D_Ebelow D2 D4])
 end
-function enlarge_southwest_corner((row, col), envs::CTMRGEnv, partfunc::InfinitePF)
-    E_south = envs.edges[SOUTH, _next(row, end), col]
-    C_southwest = envs.corners[SOUTHWEST, _next(row, end), _prev(col, end)]
-    E_west = envs.edges[WEST, row, _prev(col, end)]
+function enlarge_southwest_corner((row, col), env::CTMRGEnv, partfunc::InfinitePF)
+    E_south = env.edges[SOUTH, _next(row, end), col]
+    C_southwest = env.corners[SOUTHWEST, _next(row, end), _prev(col, end)]
+    E_west = env.edges[WEST, row, _prev(col, end)]
     return enlarge_southwest_corner(E_south, C_southwest, E_west, partfunc[row, col])
 end
 function enlarge_southwest_corner(
@@ -937,7 +937,7 @@ function renormalize_corner(
 end
 
 """
-    renormalize_northwest_corner((row, col), enlarged_envs::CTMRGEnv, P_left, P_right)
+    renormalize_northwest_corner((row, col), enlarged_env::CTMRGEnv, P_left, P_right)
     renormalize_northwest_corner(quadrant::AbstractTensorMap{T,S,3,3}, P_left, P_right) where {T,S}
     renormalize_northwest_corner(E_west, C_northwest, E_north, P_left, P_right, ket::PEPSTensor, bra::PEPSTensor=ket)
     renormalize_northwest_corner(quadrant::AbstractTensorMap{T,S,2,2}, P_left, P_right) where {T,S}
@@ -959,9 +959,9 @@ Here `A` denotes either:
 - a local pair of 'ket' and 'bra' `PEPSTensor`s
 - a `PartitionFunctionTensor`.
 """
-function renormalize_northwest_corner((row, col), enlarged_envs, P_left, P_right)
+function renormalize_northwest_corner((row, col), enlarged_env, P_left, P_right)
     return renormalize_northwest_corner(
-        enlarged_envs[NORTHWEST, row, col],
+        enlarged_env[NORTHWEST, row, col],
         P_left[NORTH, row, col],
         P_right[WEST, _next(row, end), col],
     )
@@ -1001,7 +1001,7 @@ function renormalize_northwest_corner(
 end
 
 """
-    renormalize_northeast_corner((row, col), enlarged_envs::CTMRGEnv, P_left, P_right)
+    renormalize_northeast_corner((row, col), enlarged_env::CTMRGEnv, P_left, P_right)
     renormalize_northwest_corner(quadrant::AbstractTensorMap{T,S,3,3}, P_left, P_right) where {T,S}
     renormalize_northeast_corner(E_north, C_northeast, E_east, P_left, P_right, ket::PEPSTensor, bra::PEPSTensor=ket)
     renormalize_northwest_corner(quadrant::AbstractTensorMap{T,S,2,2}, P_left, P_right) where {T,S}
@@ -1023,9 +1023,9 @@ Here `A` denotes either:
 - a local pair of 'ket' and 'bra' `PEPSTensor`s
 - a `PartitionFunctionTensor`.
 """
-function renormalize_northeast_corner((row, col), enlarged_envs, P_left, P_right)
+function renormalize_northeast_corner((row, col), enlarged_env, P_left, P_right)
     return renormalize_northeast_corner(
-        enlarged_envs[NORTHEAST, row, col],
+        enlarged_env[NORTHEAST, row, col],
         P_left[EAST, row, col],
         P_right[NORTH, row, _prev(col, end)],
     )
@@ -1067,7 +1067,7 @@ function renormalize_northeast_corner(
 end
 
 """
-    renormalize_southeast_corner((row, col), enlarged_envs::CTMRGEnv, P_left, P_right)
+    renormalize_southeast_corner((row, col), enlarged_env::CTMRGEnv, P_left, P_right)
     renormalize_southeast_corner(quadrant::AbstractTensorMap{T,S,3,3}, P_left, P_right) where {T,S}
     renormalize_southeast_corner(E_east, C_southeast, E_south, P_left, P_right, ket::PEPSTensor, bra::PEPSTensor=ket)
     renormalize_southeast_corner(quadrant::AbstractTensorMap{T,S,2,2}, P_left, P_right) where {T,S}
@@ -1089,9 +1089,9 @@ Here `A` denotes either:
 - a local pair of 'ket' and 'bra' `PEPSTensor`s
 - a `PartitionFunctionTensor`.
 """
-function renormalize_southeast_corner((row, col), enlarged_envs, P_left, P_right)
+function renormalize_southeast_corner((row, col), enlarged_env, P_left, P_right)
     return renormalize_southeast_corner(
-        enlarged_envs[SOUTHEAST, row, col],
+        enlarged_env[SOUTHEAST, row, col],
         P_left[SOUTH, row, col],
         P_right[EAST, _prev(row, end), col],
     )
@@ -1132,7 +1132,7 @@ function renormalize_southeast_corner(
 end
 
 """
-    renormalize_southwest_corner((row, col), enlarged_envs::CTMRGEnv, P_left, P_right)
+    renormalize_southwest_corner((row, col), enlarged_env::CTMRGEnv, P_left, P_right)
     renormalize_southwest_corner(quadrant::AbstractTensorMap{T,S,3,3}, P_left, P_right) where {T,S}
     renormalize_southwest_corner(E_south, C_southwest, E_west, P_left, P_right, ket::PEPSTensor, bra::PEPSTensor=ket)
     renormalize_southwest_corner(quadrant::AbstractTensorMap{T,S,2,2}, P_left, P_right) where {T,S}
@@ -1154,9 +1154,9 @@ Here `A` denotes either:
 - a pair of 'ket' and 'bra' `PEPSTensor`s
 - a `PartitionFunctionTensor`.
 """
-function renormalize_southwest_corner((row, col), enlarged_envs, P_left, P_right)
+function renormalize_southwest_corner((row, col), enlarged_env, P_left, P_right)
     return renormalize_corner(
-        enlarged_envs[SOUTHWEST, row, col],
+        enlarged_env[SOUTHWEST, row, col],
         P_left[WEST, row, col],
         P_right[SOUTH, row, _next(col, end)],
     )
@@ -1196,7 +1196,7 @@ function renormalize_southwest_corner(
 end
 
 """
-    renormalize_bottom_corner((r, c), envs, projectors)
+    renormalize_bottom_corner((r, c), env, projectors)
 
 Apply bottom projector to southwest corner and south edge.
 ```
@@ -1207,26 +1207,26 @@ Apply bottom projector to southwest corner and south edge.
 ```
 """
 function renormalize_bottom_corner(
-    (row, col), envs::CTMRGEnv{C,<:CTMRG_PEPS_EdgeTensor}, projectors
+    (row, col), env::CTMRGEnv{C,<:CTMRG_PEPS_EdgeTensor}, projectors
 ) where {C}
-    C_southwest = envs.corners[SOUTHWEST, row, _prev(col, end)]
-    E_south = envs.edges[SOUTH, row, col]
+    C_southwest = env.corners[SOUTHWEST, row, _prev(col, end)]
+    E_south = env.edges[SOUTH, row, col]
     P_bottom = projectors[1][row]
     return @autoopt @tensor corner[χ_in; χ_out] :=
         E_south[χ_in D1 D2; χ1] * C_southwest[χ1; χ2] * P_bottom[χ2 D1 D2; χ_out]
 end
 function renormalize_bottom_corner(
-    (row, col), envs::CTMRGEnv{C,<:CTMRG_PF_EdgeTensor}, projectors
+    (row, col), env::CTMRGEnv{C,<:CTMRG_PF_EdgeTensor}, projectors
 ) where {C}
-    C_southwest = envs.corners[SOUTHWEST, row, _prev(col, end)]
-    E_south = envs.edges[SOUTH, row, col]
+    C_southwest = env.corners[SOUTHWEST, row, _prev(col, end)]
+    E_south = env.edges[SOUTH, row, col]
     P_bottom = projectors[1][row]
     return @autoopt @tensor corner[χ_in; χ_out] :=
         E_south[χ_in D1; χ1] * C_southwest[χ1; χ2] * P_bottom[χ2 D1; χ_out]
 end
 
 """
-    renormalize_top_corner((row, col), envs::CTMRGEnv, projectors)
+    renormalize_top_corner((row, col), env::CTMRGEnv, projectors)
 
 Apply top projector to northwest corner and north edge.
 ```
@@ -1236,9 +1236,9 @@ Apply top projector to northwest corner and north edge.
         | 
 ```
 """
-function renormalize_top_corner((row, col), envs::CTMRGEnv, projectors)
-    C_northwest = envs.corners[NORTHWEST, row, _prev(col, end)]
-    E_north = envs.edges[NORTH, row, col]
+function renormalize_top_corner((row, col), env::CTMRGEnv, projectors)
+    C_northwest = env.corners[NORTHWEST, row, _prev(col, end)]
+    E_north = env.edges[NORTH, row, col]
     P_top = projectors[2][_next(row, end)]
     return renormalize_top_corner(C_northwest, E_north, P_top)
 end
@@ -1254,9 +1254,9 @@ end
 # edges
 
 """
-    renormalize_north_edge((row, col), envs, P_left, P_right, ket::InfinitePEPS, bra::InfinitePEPS=ket)
+    renormalize_north_edge((row, col), env, P_left, P_right, ket::InfinitePEPS, bra::InfinitePEPS=ket)
     renormalize_north_edge(E_north, P_left, P_right, ket::PEPSTensor, bra::PEPSTensor=ket)
-    renormalize_north_edge((row, col), envs, P_left, P_right, partfunc::InfinitePartitionFunction)
+    renormalize_north_edge((row, col), env, P_left, P_right, partfunc::InfinitePartitionFunction)
     renormalize_north_edge(E_north, P_left, P_right, partfunc::PartitionFunctionTensor)
 
 Absorb a local effective tensor `A` into the north edge using the given projectors and
@@ -1274,10 +1274,10 @@ Here `A` denotes either:
 - a `PartitionFunctionTensor`.
 """
 function renormalize_north_edge(
-    (row, col), envs::CTMRGEnv, P_left, P_right, ket::InfinitePEPS, bra::InfinitePEPS=ket
+    (row, col), env::CTMRGEnv, P_left, P_right, ket::InfinitePEPS, bra::InfinitePEPS=ket
 )
     return renormalize_north_edge(
-        envs.edges[NORTH, _prev(row, end), col],
+        env.edges[NORTH, _prev(row, end), col],
         P_left[NORTH, row, col],
         P_right[NORTH, row, _prev(col, end)],
         ket[row, col],
@@ -1295,10 +1295,10 @@ function renormalize_north_edge(
         P_right[χ_W; χ1 D5 D6]
 end
 function renormalize_north_edge(
-    (row, col), envs::CTMRGEnv, P_left, P_right, partfunc::InfinitePF
+    (row, col), env::CTMRGEnv, P_left, P_right, partfunc::InfinitePF
 )
     return renormalize_north_edge(
-        envs.edges[NORTH, _prev(row, end), col],
+        env.edges[NORTH, _prev(row, end), col],
         P_left[NORTH, row, col],
         P_right[NORTH, row, _prev(col, end)],
         partfunc[row, col],
@@ -1315,9 +1315,9 @@ function renormalize_north_edge(
 end
 
 """
-    renormalize_east_edge((row, col), envs, P_top, P_bottom, ket::InfinitePEPS, bra::InfinitePEPS=ket)
+    renormalize_east_edge((row, col), env, P_top, P_bottom, ket::InfinitePEPS, bra::InfinitePEPS=ket)
     renormalize_east_edge(E_east, P_top, P_bottom, ket::PEPSTensor, bra::PEPSTensor=ket)
-    renormalize_east_edge((row, col), envs, P_top, P_bottom, partfunc::InfinitePartitionFunction)
+    renormalize_east_edge((row, col), env, P_top, P_bottom, partfunc::InfinitePartitionFunction)
     renormalize_east_edge(E_east, P_top, P_bottom, partfunc::PartitionFunctionTensor)
 
 Absorb a blocal effective tensor into the east edge using the given projectors and
@@ -1338,10 +1338,10 @@ Here `A` denotes either:
 - a `PartitionFunctionTensor`.
 """
 function renormalize_east_edge(
-    (row, col), envs::CTMRGEnv, P_bottom, P_top, ket::InfinitePEPS, bra::InfinitePEPS=ket
+    (row, col), env::CTMRGEnv, P_bottom, P_top, ket::InfinitePEPS, bra::InfinitePEPS=ket
 )
     return renormalize_east_edge(
-        envs.edges[EAST, row, _next(col, end)],
+        env.edges[EAST, row, _next(col, end)],
         P_bottom[EAST, row, col, end],
         P_top[EAST, _prev(row, end), col],
         ket[row, col],
@@ -1359,10 +1359,10 @@ function renormalize_east_edge(
         P_top[χ_N; χ1 D5 D6]
 end
 function renormalize_east_edge(
-    (row, col), envs::CTMRGEnv, P_bottom, P_top, partfunc::InfinitePF
+    (row, col), env::CTMRGEnv, P_bottom, P_top, partfunc::InfinitePF
 )
     return renormalize_east_edge(
-        envs.edges[EAST, row, _next(col, end)],
+        env.edges[EAST, row, _next(col, end)],
         P_bottom[EAST, row, col, end],
         P_top[EAST, _prev(row, end), col],
         partfunc[row, col],
@@ -1379,9 +1379,9 @@ function renormalize_east_edge(
 end
 
 """
-    renormalize_south_edge((row, col), envs, P_left, P_right, ket::InfinitePEPS, bra::InfinitePEPS=ket)
+    renormalize_south_edge((row, col), env, P_left, P_right, ket::InfinitePEPS, bra::InfinitePEPS=ket)
     renormalize_south_edge(E_south, P_left, P_right, ket::PEPSTensor, bra::PEPSTensor=ket)
-    renormalize_south_edge((row, col), envs, P_left, P_right, partfunc::InfinitePartitionFunction)
+    renormalize_south_edge((row, col), env, P_left, P_right, partfunc::InfinitePartitionFunction)
     renormalize_south_edge(E_south, P_left, P_right, partfunc::PartitionFunctionTensor)
 
 Absorb a local effective tensor into the south edge using the given projectors and
@@ -1400,10 +1400,10 @@ Here `A` denotes either:
 - a `PartitionFunctionTensor`.
 """
 function renormalize_south_edge(
-    (row, col), envs::CTMRGEnv, P_left, P_right, ket::InfinitePEPS, bra::InfinitePEPS=ket
+    (row, col), env::CTMRGEnv, P_left, P_right, ket::InfinitePEPS, bra::InfinitePEPS=ket
 )
     return renormalize_south_edge(
-        envs.edges[SOUTH, _next(row, end), col],
+        env.edges[SOUTH, _next(row, end), col],
         P_left[SOUTH, row, col],
         P_right[SOUTH, row, _next(col, end)],
         ket[row, col],
@@ -1421,10 +1421,10 @@ function renormalize_south_edge(
         P_right[χ_E; χ1 D5 D6]
 end
 function renormalize_south_edge(
-    (row, col), envs::CTMRGEnv, P_left, P_right, partfunc::InfinitePF
+    (row, col), env::CTMRGEnv, P_left, P_right, partfunc::InfinitePF
 )
     return renormalize_south_edge(
-        envs.edges[SOUTH, _next(row, end), col],
+        env.edges[SOUTH, _next(row, end), col],
         P_left[SOUTH, row, col],
         P_right[SOUTH, row, _next(col, end)],
         partfunc[row, col],
@@ -1441,9 +1441,9 @@ function renormalize_south_edge(
 end
 
 """
-    renormalize_west_edge((row, col), envs, P_top, P_bottom, ket::InfinitePEPS, bra::InfinitePEPS=ket)
+    renormalize_west_edge((row, col), env, P_top, P_bottom, ket::InfinitePEPS, bra::InfinitePEPS=ket)
     renormalize_west_edge(E_west, P_top, P_bottom, ket::PEPSTensor, bra::PEPSTensor=ket)
-    renormalize_west_edge((row, col), envs, P_top, P_bottom, partfunc::InfinitePartitionFunction)
+    renormalize_west_edge((row, col), env, P_top, P_bottom, partfunc::InfinitePartitionFunction)
     renormalize_west_edge(E_west, P_top, P_bottom, partfunc::PartitionFunctionTensor)
 
 Absorb a local effective tensor into the west edge using the given projectors and
@@ -1465,14 +1465,14 @@ Here `A` denotes either:
 """
 function renormalize_west_edge(  # For simultaneous CTMRG scheme
     (row, col),
-    envs::CTMRGEnv,
+    env::CTMRGEnv,
     P_bottom::Array{Pb,3},
     P_top::Array{Pt,3},
     ket::InfinitePEPS,
     bra::InfinitePEPS=ket,
 ) where {Pt,Pb}
     return renormalize_west_edge(
-        envs.edges[WEST, row, _prev(col, end)],
+        env.edges[WEST, row, _prev(col, end)],
         P_bottom[WEST, row, col],
         P_top[WEST, _next(row, end), col],
         ket[row, col],
@@ -1481,13 +1481,13 @@ function renormalize_west_edge(  # For simultaneous CTMRG scheme
 end
 function renormalize_west_edge(  # For sequential CTMRG scheme
     (row, col),
-    envs::CTMRGEnv,
+    env::CTMRGEnv,
     projectors,
     ket::InfinitePEPS,
     bra::InfinitePEPS=ket,
 )
     return renormalize_west_edge(
-        envs.edges[WEST, row, _prev(col, end)],
+        env.edges[WEST, row, _prev(col, end)],
         projectors[1][row],
         projectors[2][_next(row, end)],
         ket[row, col],
@@ -1506,13 +1506,13 @@ function renormalize_west_edge(
 end
 function renormalize_west_edge(  # For simultaneous CTMRG scheme
     (row, col),
-    envs::CTMRGEnv,
+    env::CTMRGEnv,
     P_bottom::Array{Pb,3},
     P_top::Array{Pt,3},
     partfunc::InfinitePF,
 ) where {Pt,Pb}
     return renormalize_west_edge(
-        envs.edges[WEST, row, _prev(col, end)],
+        env.edges[WEST, row, _prev(col, end)],
         P_bottom[WEST, row, col],
         P_top[WEST, _next(row, end), col],
         partfunc[row, col],
@@ -1520,12 +1520,12 @@ function renormalize_west_edge(  # For simultaneous CTMRG scheme
 end
 function renormalize_west_edge(  # For sequential CTMRG scheme
     (row, col),
-    envs::CTMRGEnv,
+    env::CTMRGEnv,
     projectors,
     partfunc::InfinitePF,
 )
     return renormalize_west_edge(
-        envs.edges[WEST, row, _prev(col, end)],
+        env.edges[WEST, row, _prev(col, end)],
         projectors[1][row],
         projectors[2][_next(row, end)],
         partfunc[row, col],
@@ -1566,52 +1566,52 @@ function fix_gauge_corner(
 end
 
 """
-    fix_gauge_northwest_corner((row, col), envs, signs)
+    fix_gauge_northwest_corner((row, col), env, signs)
 
 Apply `fix_gauge_corner` to the northwest corner with appropriate row and column indices.
 """
-function fix_gauge_northwest_corner((row, col), envs::CTMRGEnv, signs)
+function fix_gauge_northwest_corner((row, col), env::CTMRGEnv, signs)
     return fix_gauge_corner(
-        envs.corners[NORTHWEST, row, col],
+        env.corners[NORTHWEST, row, col],
         signs[WEST, row, col],
         signs[NORTH, row, _next(col, end)],
     )
 end
 
 """
-    fix_gauge_northeast_corner((row, col), envs, signs)
+    fix_gauge_northeast_corner((row, col), env, signs)
 
 Apply `fix_gauge_corner` to the northeast corner with appropriate row and column indices.
 """
-function fix_gauge_northeast_corner((row, col), envs::CTMRGEnv, signs)
+function fix_gauge_northeast_corner((row, col), env::CTMRGEnv, signs)
     return fix_gauge_corner(
-        envs.corners[NORTHEAST, row, col],
+        env.corners[NORTHEAST, row, col],
         signs[NORTH, row, col],
         signs[EAST, _next(row, end), col],
     )
 end
 
 """
-    fix_gauge_southeast_corner((row, col), envs, signs)
+    fix_gauge_southeast_corner((row, col), env, signs)
 
 Apply `fix_gauge_corner` to the southeast corner with appropriate row and column indices.
 """
-function fix_gauge_southeast_corner((row, col), envs::CTMRGEnv, signs)
+function fix_gauge_southeast_corner((row, col), env::CTMRGEnv, signs)
     return fix_gauge_corner(
-        envs.corners[SOUTHEAST, row, col],
+        env.corners[SOUTHEAST, row, col],
         signs[EAST, row, col],
         signs[SOUTH, row, _prev(col, end)],
     )
 end
 
 """
-    fix_gauge_southwest_corner((row, col), envs, signs)
+    fix_gauge_southwest_corner((row, col), env, signs)
 
 Apply `fix_gauge_corner` to the southwest corner with appropriate row and column indices.
 """
-function fix_gauge_southwest_corner((row, col), envs::CTMRGEnv, signs)
+function fix_gauge_southwest_corner((row, col), env::CTMRGEnv, signs)
     return fix_gauge_corner(
-        envs.corners[SOUTHWEST, row, col],
+        env.corners[SOUTHWEST, row, col],
         signs[SOUTH, row, col],
         signs[WEST, _prev(row, end), col],
     )
@@ -1642,50 +1642,50 @@ function fix_gauge_edge(
 end
 
 """
-    fix_gauge_north_edge((row, col), envs, signs)
+    fix_gauge_north_edge((row, col), env, signs)
 
 Apply `fix_gauge_edge` to the north edge with appropriate row and column indices.
 """
-function fix_gauge_north_edge((row, col), envs::CTMRGEnv, signs)
+function fix_gauge_north_edge((row, col), env::CTMRGEnv, signs)
     return fix_gauge_edge(
-        envs.edges[NORTH, row, col],
+        env.edges[NORTH, row, col],
         signs[NORTH, row, col],
         signs[NORTH, row, _next(col, end)],
     )
 end
 
 """
-    fix_gauge_east_edge((row, col), envs, signs)
+    fix_gauge_east_edge((row, col), env, signs)
 
 Apply `fix_gauge_edge` to the east edge with appropriate row and column indices.
 """
-function fix_gauge_east_edge((row, col), envs::CTMRGEnv, signs)
+function fix_gauge_east_edge((row, col), env::CTMRGEnv, signs)
     return fix_gauge_edge(
-        envs.edges[EAST, row, col], signs[EAST, row, col], signs[EAST, _next(row, end), col]
+        env.edges[EAST, row, col], signs[EAST, row, col], signs[EAST, _next(row, end), col]
     )
 end
 
 """
-    fix_gauge_south_edge((row, col), envs, signs)
+    fix_gauge_south_edge((row, col), env, signs)
 
 Apply `fix_gauge_edge` to the south edge with appropriate row and column indices.
 """
-function fix_gauge_south_edge((row, col), envs::CTMRGEnv, signs)
+function fix_gauge_south_edge((row, col), env::CTMRGEnv, signs)
     return fix_gauge_edge(
-        envs.edges[SOUTH, row, col],
+        env.edges[SOUTH, row, col],
         signs[SOUTH, row, col],
         signs[SOUTH, row, _prev(col, end)],
     )
 end
 
 """
-    fix_gauge_south_edge((row, col), envs, signs)
+    fix_gauge_south_edge((row, col), env, signs)
 
 Apply `fix_gauge_edge` to the west edge with appropriate row and column indices.
 """
-function fix_gauge_west_edge((row, col), envs::CTMRGEnv, signs)
+function fix_gauge_west_edge((row, col), env::CTMRGEnv, signs)
     return fix_gauge_edge(
-        envs.edges[WEST, row, col], signs[WEST, row, col], signs[WEST, _prev(row, end), col]
+        env.edges[WEST, row, col], signs[WEST, row, col], signs[WEST, _prev(row, end), col]
     )
 end
 
diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl
index c5225f9d2..3dc5552e0 100644
--- a/src/algorithms/ctmrg/ctmrg.jl
+++ b/src/algorithms/ctmrg/ctmrg.jl
@@ -40,8 +40,9 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRGAlgorithm)
 
     return LoggingExtras.withlevel(; alg.verbosity) do
         ctmrg_loginit!(log, η, state, envinit)
+        local info
         for iter in 1:(alg.maxiter)
-            env, = ctmrg_iteration(state, env, alg)  # Grow and renormalize in all 4 directions
+            env, info = ctmrg_iteration(state, env, alg)  # Grow and renormalize in all 4 directions
             η, CS, TS = calc_convergence(env, CS, TS)
 
             if η ≤ alg.tol && iter ≥ alg.miniter
@@ -54,7 +55,7 @@ function MPSKit.leading_boundary(envinit, state, alg::CTMRGAlgorithm)
                 ctmrg_logiter!(log, iter, η, state, env)
             end
         end
-        return env
+        return env, info
     end
 end
 
@@ -100,27 +101,27 @@ function _singular_value_distance((S₁, S₂))
 end
 
 """
-    calc_convergence(envs, CS_old, TS_old)
-    calc_convergence(envs_new::CTMRGEnv, envs_old::CTMRGEnv)
+    calc_convergence(env, CS_old, TS_old)
+    calc_convergence(env_new::CTMRGEnv, env_old::CTMRGEnv)
 
-Given a new environment `envs`, compute the maximal singular value distance.
+Given a new environment `env`, compute the maximal singular value distance.
 This determined either from the previous corner and edge singular values
 `CS_old` and `TS_old`, or alternatively, directly from the old environment.
 """
-function calc_convergence(envs, CS_old, TS_old)
-    CS_new = map(x -> tsvd(x)[2], envs.corners)
+function calc_convergence(env, CS_old, TS_old)
+    CS_new = map(x -> tsvd(x)[2], env.corners)
     ΔCS = maximum(_singular_value_distance, zip(CS_old, CS_new))
 
-    TS_new = map(x -> tsvd(x)[2], envs.edges)
+    TS_new = map(x -> tsvd(x)[2], env.edges)
     ΔTS = maximum(_singular_value_distance, zip(TS_old, TS_new))
 
     @debug "maxᵢ|Cⁿ⁺¹ - Cⁿ|ᵢ = $ΔCS   maxᵢ|Tⁿ⁺¹ - Tⁿ|ᵢ = $ΔTS"
 
     return max(ΔCS, ΔTS), CS_new, TS_new
 end
-function calc_convergence(envs_new::CTMRGEnv, envs_old::CTMRGEnv)
-    CS_old = map(x -> tsvd(x)[2], envs_old.corners)
-    TS_old = map(x -> tsvd(x)[2], envs_old.edges)
-    return calc_convergence(envs_new, CS_old, TS_old)
+function calc_convergence(env_new::CTMRGEnv, env_old::CTMRGEnv)
+    CS_old = map(x -> tsvd(x)[2], env_old.corners)
+    TS_old = map(x -> tsvd(x)[2], env_old.edges)
+    return calc_convergence(env_new, CS_old, TS_old)
 end
 @non_differentiable calc_convergence(args...)
diff --git a/src/algorithms/ctmrg/projectors.jl b/src/algorithms/ctmrg/projectors.jl
index 39aa58ee7..61965f681 100644
--- a/src/algorithms/ctmrg/projectors.jl
+++ b/src/algorithms/ctmrg/projectors.jl
@@ -56,6 +56,15 @@ Projector algorithm implementing projectors from SVDing the full 4x4 CTMRG envir
     verbosity::Int = 0
 end
 
+# TODO: add `LinearAlgebra.cond` to TensorKit
+# Compute condition number smax / smin for diagonal singular value TensorMap
+function _condition_number(S::AbstractTensorMap)
+    smax = maximum(first ∘ last, blocks(S))
+    smin = maximum(last ∘ last, blocks(S))
+    return smax / smin
+end
+@non_differentiable _condition_number(S::AbstractTensorMap)
+
 """
     compute_projector(enlarged_corners, coordinate, alg::ProjectorAlgorithm)
 
@@ -66,16 +75,20 @@ function compute_projector(enlarged_corners, coordinate, alg::HalfInfiniteProjec
     # SVD half-infinite environment
     halfinf = half_infinite_environment(enlarged_corners...)
     svd_alg = svd_algorithm(alg, coordinate)
-    U, S, V, err = PEPSKit.tsvd!(halfinf, svd_alg; trunc=alg.trscheme)
-    # Compute SVD truncation error and check for degenerate singular values
+    U, S, V, truncation_error = PEPSKit.tsvd!(halfinf, svd_alg; trunc=alg.trscheme)
+
+    # Check for degenerate singular values
     Zygote.isderiving() && ignore_derivatives() do
         if alg.verbosity > 0 && is_degenerate_spectrum(S)
             svals = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(S))
             @warn("degenerate singular values detected: ", svals)
         end
     end
+
     P_left, P_right = contract_projectors(U, S, V, enlarged_corners...)
-    return (P_left, P_right), (; err, U, S, V)
+    truncation_error /= norm(S)
+    condition_number = @ignore_derivatives(_condition_number(S))
+    return (P_left, P_right), (; truncation_error, condition_number, U, S, V)
 end
 function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjector)
     halfinf_left = half_infinite_environment(enlarged_corners[1], enlarged_corners[2])
@@ -84,14 +97,17 @@ function compute_projector(enlarged_corners, coordinate, alg::FullInfiniteProjec
     # SVD full-infinite environment
     fullinf = full_infinite_environment(halfinf_left, halfinf_right)
     svd_alg = svd_algorithm(alg, coordinate)
-    U, S, V, err = PEPSKit.tsvd!(fullinf, svd_alg; trunc=alg.trscheme)
-    # Compute SVD truncation error and check for degenerate singular values
+    U, S, V, truncation_error = PEPSKit.tsvd!(fullinf, svd_alg; trunc=alg.trscheme)
+
+    # Check for degenerate singular values
     Zygote.isderiving() && ignore_derivatives() do
         if alg.verbosity > 0 && is_degenerate_spectrum(S)
             svals = TensorKit.SectorDict(c => diag(b) for (c, b) in blocks(S))
             @warn("degenerate singular values detected: ", svals)
         end
     end
+
     P_left, P_right = contract_projectors(U, S, V, halfinf_left, halfinf_right)
-    return (P_left, P_right), (; err, U, S, V)
+    condition_number = @ignore_derivatives(_condition_number(S))
+    return (P_left, P_right), (; truncation_error, condition_number, U, S, V)
 end
diff --git a/src/algorithms/ctmrg/sequential.jl b/src/algorithms/ctmrg/sequential.jl
index 4404a7f5a..4d50a960f 100644
--- a/src/algorithms/ctmrg/sequential.jl
+++ b/src/algorithms/ctmrg/sequential.jl
@@ -30,103 +30,101 @@ function SequentialCTMRG(;
     )
 end
 
-function ctmrg_iteration(state, envs::CTMRGEnv, alg::SequentialCTMRG)
-    ϵ = zero(real(scalartype(state)))
+function ctmrg_iteration(state, env::CTMRGEnv, alg::SequentialCTMRG)
+    truncation_error = zero(real(scalartype(state)))
+    condition_number = zero(real(scalartype(state)))
     for _ in 1:4 # rotate
         for col in 1:size(state, 2) # left move column-wise
-            projectors, info = sequential_projectors(col, state, envs, alg.projector_alg)
-            envs = renormalize_sequentially(col, projectors, state, envs)
-            ϵ = max(ϵ, info.err)
+            projectors, info = sequential_projectors(col, state, env, alg.projector_alg)
+            env = renormalize_sequentially(col, projectors, state, env)
+            truncation_error = max(truncation_error, info.truncation_error)
+            condition_number = max(condition_number, info.condition_number)
         end
         state = rotate_north(state, EAST)
-        envs = rotate_north(envs, EAST)
+        env = rotate_north(env, EAST)
     end
 
-    return envs, (; err=ϵ)
+    return env, (; truncation_error, condition_number)
 end
 
 """
-    sequential_projectors(col::Int, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgorithm)
-    sequential_projectors(coordinate::NTuple{3,Int}, state::InfinitePEPS, envs::CTMRGEnv, alg::ProjectorAlgorithm)
+    sequential_projectors(col::Int, state::InfinitePEPS, env::CTMRGEnv, alg::ProjectorAlgorithm)
+    sequential_projectors(coordinate::NTuple{3,Int}, state::InfinitePEPS, env::CTMRGEnv, alg::ProjectorAlgorithm)
 
 Compute CTMRG projectors in the `:sequential` scheme either for an entire column `col` or
 for a specific `coordinate` (where `dir=WEST` is already implied in the `:sequential` scheme).
 """
 function sequential_projectors(
-    col::Int, state::InfiniteSquareNetwork, envs::CTMRGEnv, alg::ProjectorAlgorithm
+    col::Int, state::InfiniteSquareNetwork, env::CTMRGEnv, alg::ProjectorAlgorithm
 )
-    # SVD half-infinite environment column-wise
-    ϵ = Zygote.Buffer(zeros(real(scalartype(envs)), size(envs, 2)))
-    coordinates = eachcoordinate(envs)[:, col]
-    projectors = dtmap(coordinates) do (r, c)
-        trscheme = truncation_scheme(alg, envs.edges[WEST, _prev(r, size(envs, 2)), c])
+    coordinates = eachcoordinate(env)[:, col]
+    proj_and_info = dtmap(coordinates) do (r, c)
+        trscheme = truncation_scheme(alg, env.edges[WEST, _prev(r, size(env, 2)), c])
         proj, info = sequential_projectors(
-            (WEST, r, c), state, envs, @set(alg.trscheme = trscheme)
+            (WEST, r, c), state, env, @set(alg.trscheme = trscheme)
         )
-        ϵ[r] = info.err / norm(info.S)
-        return proj
+        return proj, info
     end
-
-    return (map(first, projectors), map(last, projectors)), (; err=maximum(copy(ϵ)))
+    return _split_proj_and_info(proj_and_info)
 end
 function sequential_projectors(
     coordinate::NTuple{3,Int},
     state::InfiniteSquareNetwork,
-    envs::CTMRGEnv,
+    env::CTMRGEnv,
     alg::HalfInfiniteProjector,
 )
     _, r, c = coordinate
-    r′ = _prev(r, size(envs, 2))
-    Q1 = TensorMap(EnlargedCorner(state, envs, (SOUTHWEST, r, c)), SOUTHWEST)
-    Q2 = TensorMap(EnlargedCorner(state, envs, (NORTHWEST, r′, c)), NORTHWEST)
+    r′ = _prev(r, size(env, 2))
+    Q1 = TensorMap(EnlargedCorner(state, env, (SOUTHWEST, r, c)), SOUTHWEST)
+    Q2 = TensorMap(EnlargedCorner(state, env, (NORTHWEST, r′, c)), NORTHWEST)
     return compute_projector((Q1, Q2), coordinate, alg)
 end
 function sequential_projectors(
     coordinate::NTuple{3,Int},
     state::InfiniteSquareNetwork,
-    envs::CTMRGEnv,
+    env::CTMRGEnv,
     alg::FullInfiniteProjector,
 )
-    rowsize, colsize = size(envs)[2:3]
+    rowsize, colsize = size(env)[2:3]
     coordinate_nw = _next_coordinate(coordinate, rowsize, colsize)
     coordinate_ne = _next_coordinate(coordinate_nw, rowsize, colsize)
     coordinate_se = _next_coordinate(coordinate_ne, rowsize, colsize)
     ec = (
-        TensorMap(EnlargedCorner(state, envs, coordinate_se), SOUTHEAST),
-        TensorMap(EnlargedCorner(state, envs, coordinate), SOUTHWEST),
-        TensorMap(EnlargedCorner(state, envs, coordinate_nw), NORTHWEST),
-        TensorMap(EnlargedCorner(state, envs, coordinate_ne), NORTHEAST),
+        TensorMap(EnlargedCorner(state, env, coordinate_se), SOUTHEAST),
+        TensorMap(EnlargedCorner(state, env, coordinate), SOUTHWEST),
+        TensorMap(EnlargedCorner(state, env, coordinate_nw), NORTHWEST),
+        TensorMap(EnlargedCorner(state, env, coordinate_ne), NORTHEAST),
     )
     return compute_projector(ec, coordinate, alg)
 end
 
 """
-    renormalize_sequentially(col::Int, projectors, state, envs)
+    renormalize_sequentially(col::Int, projectors, state, env)
 
 Renormalize one column of the CTMRG environment.
 """
-function renormalize_sequentially(col::Int, projectors, state, envs)
-    corners = Zygote.Buffer(envs.corners)
-    edges = Zygote.Buffer(envs.edges)
+function renormalize_sequentially(col::Int, projectors, state, env)
+    corners = Zygote.Buffer(env.corners)
+    edges = Zygote.Buffer(env.edges)
 
     for (dir, r, c) in eachcoordinate(state, 1:4)
         (c == col && dir in [SOUTHWEST, NORTHWEST]) && continue
-        corners[dir, r, c] = envs.corners[dir, r, c]
+        corners[dir, r, c] = env.corners[dir, r, c]
     end
     for (dir, r, c) in eachcoordinate(state, 1:4)
         (c == col && dir == WEST) && continue
-        edges[dir, r, c] = envs.edges[dir, r, c]
+        edges[dir, r, c] = env.edges[dir, r, c]
     end
 
     # Apply projectors to renormalize corners and edge
-    for row in axes(envs.corners, 2)
-        C_southwest = renormalize_bottom_corner((row, col), envs, projectors)
+    for row in axes(env.corners, 2)
+        C_southwest = renormalize_bottom_corner((row, col), env, projectors)
         corners[SOUTHWEST, row, col] = C_southwest / norm(C_southwest)
 
-        C_northwest = renormalize_top_corner((row, col), envs, projectors)
+        C_northwest = renormalize_top_corner((row, col), env, projectors)
         corners[NORTHWEST, row, col] = C_northwest / norm(C_northwest)
 
-        E_west = renormalize_west_edge((row, col), envs, projectors, state)
+        E_west = renormalize_west_edge((row, col), env, projectors, state)
         edges[WEST, row, col] = E_west / norm(E_west)
     end
 
diff --git a/src/algorithms/ctmrg/simultaneous.jl b/src/algorithms/ctmrg/simultaneous.jl
index c592d64e2..79b8b6e1e 100644
--- a/src/algorithms/ctmrg/simultaneous.jl
+++ b/src/algorithms/ctmrg/simultaneous.jl
@@ -30,53 +30,47 @@ function SimultaneousCTMRG(;
     )
 end
 
-function ctmrg_iteration(state, envs::CTMRGEnv, alg::SimultaneousCTMRG)
+function ctmrg_iteration(state, env::CTMRGEnv, alg::SimultaneousCTMRG)
     enlarged_corners = dtmap(eachcoordinate(state, 1:4)) do idx
-        return TensorMap(EnlargedCorner(state, envs, idx), idx[1])
+        return TensorMap(EnlargedCorner(state, env, idx), idx[1])
     end  # expand environment
-    projectors, info = simultaneous_projectors(enlarged_corners, envs, alg.projector_alg)  # compute projectors on all coordinates
-    envs′ = renormalize_simultaneously(enlarged_corners, projectors, state, envs)  # renormalize enlarged corners
-    return envs′, info
+    projectors, info = simultaneous_projectors(enlarged_corners, env, alg.projector_alg)  # compute projectors on all coordinates
+    env′ = renormalize_simultaneously(enlarged_corners, projectors, state, env)  # renormalize enlarged corners
+    return env′, info
 end
 
-# Pre-allocate U, S, and V tensor as Zygote buffers to make it differentiable
-function _prealloc_svd(edges::Array{E,N}) where {E,N}
-    Sc = scalartype(E)
-    U = Zygote.Buffer(map(e -> zeros(Sc, space(e)), edges))
-    V = Zygote.Buffer(map(e -> zeros(Sc, domain(e), codomain(e)), edges))
-    S = Zygote.Buffer(U.data, tensormaptype(spacetype(E), 1, 1, real(Sc)))  # Corner type but with real numbers
-    return U, S, V
+# Work-around to stop Zygote from choking on first execution (sometimes)
+# Split up map returning projectors and info into separate arrays
+function _split_proj_and_info(proj_and_info)
+    P_left = map(x -> x[1][1], proj_and_info)
+    P_right = map(x -> x[1][2], proj_and_info)
+    truncation_error = maximum(x -> x[2].truncation_error, proj_and_info)
+    condition_number = maximum(x -> x[2].condition_number, proj_and_info)
+    U = map(x -> x[2].U, proj_and_info)
+    S = map(x -> x[2].S, proj_and_info)
+    V = map(x -> x[2].V, proj_and_info)
+    info = (; truncation_error, condition_number, U, S, V)
+    return (P_left, P_right), info
 end
 
 """
-    simultaneous_projectors(enlarged_corners::Array{E,3}, envs::CTMRGEnv, alg::ProjectorAlgorithm)
+    simultaneous_projectors(enlarged_corners::Array{E,3}, env::CTMRGEnv, alg::ProjectorAlgorithm)
     simultaneous_projectors(coordinate, enlarged_corners::Array{E,3}, alg::ProjectorAlgorithm)
 
 Compute CTMRG projectors in the `:simultaneous` scheme either for all provided
 enlarged corners or on a specific `coordinate`.
 """
 function simultaneous_projectors(
-    enlarged_corners::Array{E,3}, envs::CTMRGEnv, alg::ProjectorAlgorithm
+    enlarged_corners::Array{E,3}, env::CTMRGEnv, alg::ProjectorAlgorithm
 ) where {E}
-    U, S, V = _prealloc_svd(envs.edges)
-    ϵ = Zygote.Buffer(zeros(real(scalartype(envs)), size(envs)))
-
-    projectors = dtmap(eachcoordinate(envs, 1:4)) do coordinate
-        coordinate′ = _next_coordinate(coordinate, size(envs)[2:3]...)
-        trscheme = truncation_scheme(alg, envs.edges[coordinate[1], coordinate′[2:3]...])
-        proj, info = simultaneous_projectors(
+    proj_and_info = dtmap(eachcoordinate(env, 1:4)) do coordinate
+        coordinate′ = _next_coordinate(coordinate, size(env)[2:3]...)
+        trscheme = truncation_scheme(alg, env.edges[coordinate[1], coordinate′[2:3]...])
+        return simultaneous_projectors(
             coordinate, enlarged_corners, @set(alg.trscheme = trscheme)
         )
-        U[coordinate...] = info.U
-        S[coordinate...] = info.S
-        V[coordinate...] = info.V
-        ϵ[coordinate...] = info.err / norm(info.S)
-        return proj
     end
-
-    P_left = map(first, projectors)
-    P_right = map(last, projectors)
-    return (P_left, P_right), (; err=maximum(copy(ϵ)), U=copy(U), S=copy(S), V=copy(V))
+    return _split_proj_and_info(proj_and_info)
 end
 function simultaneous_projectors(
     coordinate, enlarged_corners::Array{E,3}, alg::HalfInfiniteProjector
@@ -102,26 +96,26 @@ function simultaneous_projectors(
 end
 
 """
-    renormalize_simultaneously(enlarged_corners, projectors, state, envs)
+    renormalize_simultaneously(enlarged_corners, projectors, state, env)
 
 Renormalize all enlarged corners and edges simultaneously.
 """
-function renormalize_simultaneously(enlarged_corners, projectors, state, envs)
+function renormalize_simultaneously(enlarged_corners, projectors, state, env)
     P_left, P_right = projectors
-    coordinates = eachcoordinate(envs, 1:4)
+    coordinates = eachcoordinate(env, 1:4)
     corners_edges = dtmap(coordinates) do (dir, r, c)
         if dir == NORTH
             corner = renormalize_northwest_corner((r, c), enlarged_corners, P_left, P_right)
-            edge = renormalize_north_edge((r, c), envs, P_left, P_right, state)
+            edge = renormalize_north_edge((r, c), env, P_left, P_right, state)
         elseif dir == EAST
             corner = renormalize_northeast_corner((r, c), enlarged_corners, P_left, P_right)
-            edge = renormalize_east_edge((r, c), envs, P_left, P_right, state)
+            edge = renormalize_east_edge((r, c), env, P_left, P_right, state)
         elseif dir == SOUTH
             corner = renormalize_southeast_corner((r, c), enlarged_corners, P_left, P_right)
-            edge = renormalize_south_edge((r, c), envs, P_left, P_right, state)
+            edge = renormalize_south_edge((r, c), env, P_left, P_right, state)
         elseif dir == WEST
             corner = renormalize_southwest_corner((r, c), enlarged_corners, P_left, P_right)
-            edge = renormalize_west_edge((r, c), envs, P_left, P_right, state)
+            edge = renormalize_west_edge((r, c), env, P_left, P_right, state)
         end
         return corner / norm(corner), edge / norm(edge)
     end
diff --git a/src/algorithms/ctmrg/sparse_environments.jl b/src/algorithms/ctmrg/sparse_environments.jl
index cbd1a660a..616e89f1e 100644
--- a/src/algorithms/ctmrg/sparse_environments.jl
+++ b/src/algorithms/ctmrg/sparse_environments.jl
@@ -29,39 +29,39 @@ function EnlargedCorner(
 end
 
 """
-    EnlargedCorner(network::InfiniteSquareNetwork, envs, coordinates)
+    EnlargedCorner(network::InfiniteSquareNetwork, env, coordinates)
 
 Construct an enlarged corner with the correct row and column indices based on the given
 `coordinates` which are of the form `(dir, row, col)`.
 """
-function EnlargedCorner(network::InfiniteSquareNetwork, envs, coordinates)
+function EnlargedCorner(network::InfiniteSquareNetwork, env, coordinates)
     dir, r, c = coordinates
     if dir == NORTHWEST
         return EnlargedCorner(
-            envs.corners[NORTHWEST, _prev(r, end), _prev(c, end)],
-            envs.edges[WEST, r, _prev(c, end)],
-            envs.edges[NORTH, _prev(r, end), c],
+            env.corners[NORTHWEST, _prev(r, end), _prev(c, end)],
+            env.edges[WEST, r, _prev(c, end)],
+            env.edges[NORTH, _prev(r, end), c],
             network[r, c],
         )
     elseif dir == NORTHEAST
         return EnlargedCorner(
-            envs.corners[NORTHEAST, _prev(r, end), _next(c, end)],
-            envs.edges[NORTH, _prev(r, end), c],
-            envs.edges[EAST, r, _next(c, end)],
+            env.corners[NORTHEAST, _prev(r, end), _next(c, end)],
+            env.edges[NORTH, _prev(r, end), c],
+            env.edges[EAST, r, _next(c, end)],
             network[r, c],
         )
     elseif dir == SOUTHEAST
         return EnlargedCorner(
-            envs.corners[SOUTHEAST, _next(r, end), _next(c, end)],
-            envs.edges[EAST, r, _next(c, end)],
-            envs.edges[SOUTH, _next(r, end), c],
+            env.corners[SOUTHEAST, _next(r, end), _next(c, end)],
+            env.edges[EAST, r, _next(c, end)],
+            env.edges[SOUTH, _next(r, end), c],
             network[r, c],
         )
     elseif dir == SOUTHWEST
         return EnlargedCorner(
-            envs.corners[SOUTHWEST, _next(r, end), _prev(c, end)],
-            envs.edges[SOUTH, _next(r, end), c],
-            envs.edges[WEST, r, _prev(c, end)],
+            env.corners[SOUTHWEST, _next(r, end), _prev(c, end)],
+            env.edges[SOUTH, _next(r, end), c],
+            env.edges[WEST, r, _prev(c, end)],
             network[r, c],
         )
     end
diff --git a/src/algorithms/peps_opt.jl b/src/algorithms/optimization/fixed_point_differentiation.jl
similarity index 55%
rename from src/algorithms/peps_opt.jl
rename to src/algorithms/optimization/fixed_point_differentiation.jl
index 9d7b20568..4194143a5 100644
--- a/src/algorithms/peps_opt.jl
+++ b/src/algorithms/optimization/fixed_point_differentiation.jl
@@ -76,137 +76,6 @@ function LinSolver(;
     return LinSolver{iterscheme}(solver)
 end
 
-"""
-    PEPSOptimize{G}(; boundary_alg=Defaults.ctmrg_alg, optimizer::OptimKit.OptimizationAlgorithm=Defaults.optimizer
-                    reuse_env::Bool=true, gradient_alg::G=Defaults.gradient_alg)
-
-Algorithm struct that represent PEPS ground-state optimization using AD.
-Set the algorithm to contract the infinite PEPS in `boundary_alg`;
-currently only `CTMRGAlgorithm`s are supported. The `optimizer` computes the gradient directions
-based on the CTMRG gradient and updates the PEPS parameters. In this optimization,
-the CTMRG runs can be started on the converged environments of the previous optimizer
-step by setting `reuse_env` to true. Otherwise a random environment is used at each
-step. The CTMRG gradient itself is computed using the `gradient_alg` algorithm.
-"""
-struct PEPSOptimize{G}
-    boundary_alg::CTMRGAlgorithm
-    optimizer::OptimKit.OptimizationAlgorithm
-    reuse_env::Bool
-    gradient_alg::G
-
-    function PEPSOptimize(  # Inner constructor to prohibit illegal setting combinations
-        boundary_alg::CTMRGAlgorithm,
-        optimizer,
-        reuse_env,
-        gradient_alg::G,
-    ) where {G}
-        if gradient_alg isa GradMode
-            if boundary_alg isa SequentialCTMRG && iterscheme(gradient_alg) === :fixed
-                throw(ArgumentError(":sequential and :fixed are not compatible"))
-            end
-        end
-        return new{G}(boundary_alg, optimizer, reuse_env, gradient_alg)
-    end
-end
-function PEPSOptimize(;
-    boundary_alg=Defaults.ctmrg_alg,
-    optimizer=Defaults.optimizer,
-    reuse_env=Defaults.reuse_env,
-    gradient_alg=Defaults.gradient_alg,
-)
-    return PEPSOptimize(boundary_alg, optimizer, reuse_env, gradient_alg)
-end
-
-"""
-    fixedpoint(ψ₀::InfinitePEPS{T}, H, alg::PEPSOptimize, [env₀::CTMRGEnv];
-               finalize!=OptimKit._finalize!, symmetrization=nothing) where {T}
-    
-Optimize `ψ₀` with respect to the Hamiltonian `H` according to the parameters supplied
-in `alg`. The initial environment `env₀` serves as an initial guess for the first CTMRG run.
-By default, a random initial environment is used.
-
-The `finalize!` kwarg can be used to insert a function call after each optimization step
-by utilizing the `finalize!` kwarg of `OptimKit.optimize`.
-The function maps `(peps, envs), f, g = finalize!((peps, envs), f, g, numiter)`.
-The `symmetrization` kwarg accepts `nothing` or a `SymmetrizationStyle`, in which case the
-PEPS and PEPS gradient are symmetrized after each optimization iteration. Note that this
-requires a symmmetric `ψ₀` and `env₀` to converge properly.
-
-The function returns a `NamedTuple` which contains the following entries:
-- `peps`: final `InfinitePEPS`
-- `env`: `CTMRGEnv` corresponding to the final PEPS
-- `E`: final energy
-- `E_history`: convergence history of the energy function
-- `grad`: final energy gradient
-- `gradnorm_history`: convergence history of the energy gradient norms
-- `numfg`: total number of calls to the energy function
-"""
-function fixedpoint(
-    ψ₀::InfinitePEPS{F},
-    H,
-    alg::PEPSOptimize,
-    env₀::CTMRGEnv=CTMRGEnv(ψ₀, field(F)^20);
-    (finalize!)=OptimKit._finalize!,
-    symmetrization=nothing,
-) where {F}
-    if isnothing(symmetrization)
-        retract = peps_retract
-    else
-        retract, symm_finalize! = symmetrize_retract_and_finalize!(symmetrization)
-        fin! = finalize!  # Previous finalize!
-        finalize! = (x, f, g, numiter) -> fin!(symm_finalize!(x, f, g, numiter)..., numiter)
-    end
-
-    if scalartype(env₀) <: Real
-        env₀ = complex(env₀)
-        @warn "the provided real environment was converted to a complex environment since \
-        :fixed mode generally produces complex gauges; use :diffgauge mode instead to work \
-        with purely real environments"
-    end
-
-    (peps, env), E, ∂E, numfg, convhistory = optimize(
-        (ψ₀, env₀), alg.optimizer; retract, inner=real_inner, finalize!
-    ) do (peps, envs)
-        E, gs = withgradient(peps) do ψ
-            envs´ = hook_pullback(
-                leading_boundary,
-                envs,
-                ψ,
-                alg.boundary_alg;
-                alg_rrule=alg.gradient_alg,
-            )
-            ignore_derivatives() do
-                alg.reuse_env && update!(envs, envs´)
-            end
-            return costfun(ψ, envs´, H)
-        end
-        g = only(gs)  # `withgradient` returns tuple of gradients `gs`
-        return E, g
-    end
-
-    return (;
-        peps,
-        env,
-        E,
-        E_history=convhistory[:, 1],
-        grad=∂E,
-        gradnorm_history=convhistory[:, 2],
-        numfg,
-    )
-end
-
-# Update PEPS unit cell in non-mutating way
-# Note: Both x and η are InfinitePEPS during optimization
-function peps_retract(x, η, α)
-    peps = deepcopy(x[1])
-    peps.A .+= η.A .* α
-    env = deepcopy(x[2])
-    return (peps, env), η
-end
-
-# Take real valued part of dot product
-real_inner(_, η₁, η₂) = real(dot(η₁, η₂))
-
 #=
 Evaluating the gradient of the cost function for CTMRG:
 - The gradient of the cost function for CTMRG can be computed using automatic differentiation (AD) or explicit evaluation of the geometric sum.
@@ -222,24 +91,24 @@ function _rrule(
     state,
     alg::CTMRGAlgorithm,
 )
-    envs = leading_boundary(envinit, state, alg)
+    env, info = leading_boundary(envinit, state, alg)
 
-    function leading_boundary_diffgauge_pullback(Δenvs′)
-        Δenvs = unthunk(Δenvs′)
+    function leading_boundary_diffgauge_pullback((Δenv′, Δinfo))
+        Δenv = unthunk(Δenv′)
 
         # find partial gradients of gauge_fixed single CTMRG iteration
         f(A, x) = gauge_fix(x, ctmrg_iteration(A, x, alg)[1])[1]
-        _, env_vjp = rrule_via_ad(config, f, state, envs)
+        _, env_vjp = rrule_via_ad(config, f, state, env)
 
         # evaluate the geometric sum
         ∂f∂A(x)::typeof(state) = env_vjp(x)[2]
-        ∂f∂x(x)::typeof(envs) = env_vjp(x)[3]
-        ∂F∂envs = fpgrad(Δenvs, ∂f∂x, ∂f∂A, Δenvs, gradmode)
+        ∂f∂x(x)::typeof(env) = env_vjp(x)[3]
+        ∂F∂env = fpgrad(Δenv, ∂f∂x, ∂f∂A, Δenv, gradmode)
 
-        return NoTangent(), ZeroTangent(), ∂F∂envs, NoTangent()
+        return NoTangent(), ZeroTangent(), ∂F∂env, NoTangent()
     end
 
-    return envs, leading_boundary_diffgauge_pullback
+    return (env, info), leading_boundary_diffgauge_pullback
 end
 
 # Here f is differentiated from an pre-computed SVD with fixed U, S and V
@@ -252,9 +121,9 @@ function _rrule(
     alg::SimultaneousCTMRG,
 )
     @assert !isnothing(alg.projector_alg.svd_alg.rrule_alg)
-    envs = leading_boundary(envinit, state, alg)
-    envs_conv, info = ctmrg_iteration(state, envs, alg)
-    envs_fixed, signs = gauge_fix(envs, envs_conv)
+    env, = leading_boundary(envinit, state, alg)
+    env_conv, info = ctmrg_iteration(state, env, alg)
+    env_fixed, signs = gauge_fix(env, env_conv)
 
     # Fix SVD
     U_fixed, V_fixed = fix_relative_phases(info.U, info.V, signs)
@@ -265,21 +134,21 @@ function _rrule(
     alg_fixed = @set alg.projector_alg.svd_alg = svd_alg_fixed
     alg_fixed = @set alg_fixed.projector_alg.trscheme = notrunc()
 
-    function leading_boundary_fixed_pullback(Δenvs′)
-        Δenvs = unthunk(Δenvs′)
+    function leading_boundary_fixed_pullback((Δenv′, Δinfo))
+        Δenv = unthunk(Δenv′)
 
         f(A, x) = fix_global_phases(x, ctmrg_iteration(A, x, alg_fixed)[1])
-        _, env_vjp = rrule_via_ad(config, f, state, envs_fixed)
+        _, env_vjp = rrule_via_ad(config, f, state, env_fixed)
 
         # evaluate the geometric sum
         ∂f∂A(x)::typeof(state) = env_vjp(x)[2]
-        ∂f∂x(x)::typeof(envs) = env_vjp(x)[3]
-        ∂F∂envs = fpgrad(Δenvs, ∂f∂x, ∂f∂A, Δenvs, gradmode)
+        ∂f∂x(x)::typeof(env) = env_vjp(x)[3]
+        ∂F∂env = fpgrad(Δenv, ∂f∂x, ∂f∂A, Δenv, gradmode)
 
-        return NoTangent(), ZeroTangent(), ∂F∂envs, NoTangent()
+        return NoTangent(), ZeroTangent(), ∂F∂env, NoTangent()
     end
 
-    return envs_fixed, leading_boundary_fixed_pullback
+    return (env_fixed, info), leading_boundary_fixed_pullback
 end
 
 @doc """
diff --git a/src/algorithms/optimization/peps_optimization.jl b/src/algorithms/optimization/peps_optimization.jl
new file mode 100644
index 000000000..aecb2132b
--- /dev/null
+++ b/src/algorithms/optimization/peps_optimization.jl
@@ -0,0 +1,177 @@
+"""
+    PEPSOptimize{G}(; boundary_alg=Defaults.ctmrg_alg, optimizer::OptimKit.OptimizationAlgorithm=Defaults.optimizer
+                    reuse_env::Bool=true, gradient_alg::G=Defaults.gradient_alg)
+
+Algorithm struct that represent PEPS ground-state optimization using AD.
+Set the algorithm to contract the infinite PEPS in `boundary_alg`;
+currently only `CTMRGAlgorithm`s are supported. The `optimizer` computes the gradient directions
+based on the CTMRG gradient and updates the PEPS parameters. In this optimization,
+the CTMRG runs can be started on the converged environments of the previous optimizer
+step by setting `reuse_env` to true. Otherwise a random environment is used at each
+step. The CTMRG gradient itself is computed using the `gradient_alg` algorithm.
+"""
+struct PEPSOptimize{G}
+    boundary_alg::CTMRGAlgorithm
+    gradient_alg::G
+    optimizer::OptimKit.OptimizationAlgorithm
+    reuse_env::Bool
+    symmetrization::Union{Nothing,SymmetrizationStyle}
+
+    function PEPSOptimize(  # Inner constructor to prohibit illegal setting combinations
+        boundary_alg::CTMRGAlgorithm,
+        gradient_alg::G,
+        optimizer,
+        reuse_env,
+        symmetrization,
+    ) where {G}
+        if gradient_alg isa GradMode
+            if boundary_alg isa SequentialCTMRG && iterscheme(gradient_alg) === :fixed
+                throw(ArgumentError(":sequential and :fixed are not compatible"))
+            end
+        end
+        return new{G}(boundary_alg, gradient_alg, optimizer, reuse_env, symmetrization)
+    end
+end
+function PEPSOptimize(;
+    boundary_alg=Defaults.ctmrg_alg,
+    gradient_alg=Defaults.gradient_alg,
+    optimizer=Defaults.optimizer,
+    reuse_env=Defaults.reuse_env,
+    symmetrization=nothing,
+)
+    return PEPSOptimize(boundary_alg, gradient_alg, optimizer, reuse_env, symmetrization)
+end
+
+"""
+    fixedpoint(operator, peps₀::InfinitePEPS, env₀::CTMRGEnv; kwargs...)
+    fixedpoint(operator, peps₀::InfinitePEPS, env₀::CTMRGEnv, alg::PEPSOptimize;
+               finalize!=OptimKit._finalize!)
+    
+Find the fixed point of `operator` (i.e. the ground state) starting from `peps₀` according
+to the optimization parameters supplied in `alg`. The initial environment `env₀` serves as
+an initial guess for the first CTMRG run. By default, a random initial environment is used.
+
+The `finalize!` kwarg can be used to insert a function call after each optimization step
+by utilizing the `finalize!` kwarg of `OptimKit.optimize`.
+The function maps `(peps, env), f, g = finalize!((peps, env), f, g, numiter)`.
+The `symmetrization` kwarg accepts `nothing` or a `SymmetrizationStyle`, in which case the
+PEPS and PEPS gradient are symmetrized after each optimization iteration. Note that this
+requires a symmmetric `peps₀` and `env₀` to converge properly.
+
+The function returns the final PEPS, CTMRG environment and cost value, as well as an
+information `NamedTuple` which contains the following entries:
+- `last_gradient`: last gradient of the cost function
+- `fg_evaluations`: number of evaluations of the cost and gradient function
+- `costs`: history of cost values
+- `gradnorms`: history of gradient norms
+- `truncation_errors`: history of truncation errors of the boundary algorithm
+- `condition_numbers`: history of condition numbers of the CTMRG environments
+- `gradnorms_unitcell`: history of gradient norms for each respective unit cell entry
+- `times`: history of times each optimization step took
+"""
+function fixedpoint(operator, peps₀::InfinitePEPS, env₀::CTMRGEnv; kwargs...)
+    throw(error("method not yet implemented"))
+    alg = fixedpoint_selector(; kwargs...) # TODO: implement fixedpoint_selector
+    return fixedpoint(operator, peps₀, env₀, alg)
+end
+function fixedpoint(
+    operator,
+    peps₀::InfinitePEPS,
+    env₀::CTMRGEnv,
+    alg::PEPSOptimize;
+    (finalize!)=OptimKit._finalize!,
+)
+    # setup retract and finalize! for symmetrization
+    if isnothing(alg.symmetrization)
+        retract = peps_retract
+    else
+        retract, symm_finalize! = symmetrize_retract_and_finalize!(alg.symmetrization)
+        fin! = finalize!  # Previous finalize!
+        finalize! = (x, f, g, numiter) -> fin!(symm_finalize!(x, f, g, numiter)..., numiter)
+    end
+
+    # check realness compatibility
+    if scalartype(env₀) <: Real && iterscheme(alg.gradient_alg) == :fixed
+        env₀ = complex(env₀)
+        @warn "the provided real environment was converted to a complex environment since \
+        :fixed mode generally produces complex gauges; use :diffgauge mode instead to work \
+        with purely real environments"
+    end
+
+    # initialize info collection vectors
+    T = promote_type(real(scalartype(peps₀)), real(scalartype(env₀)))
+    truncation_errors = Vector{T}()
+    condition_numbers = Vector{T}()
+    gradnorms_unitcell = Vector{Matrix{T}}()
+    times = Vector{Float64}()
+
+    # optimize operator cost function
+    (peps_final, env_final), cost, ∂cost, numfg, convergence_history = optimize(
+        (peps₀, env₀), alg.optimizer; retract, inner=real_inner, finalize!
+    ) do (peps, env)
+        start_time = time_ns()
+        E, gs = withgradient(peps) do ψ
+            env′, info = hook_pullback(
+                leading_boundary,
+                env,
+                ψ,
+                alg.boundary_alg;
+                alg_rrule=alg.gradient_alg,
+            )
+            ignore_derivatives() do
+                alg.reuse_env && update!(env, env′)
+                push!(truncation_errors, info.truncation_error)
+                push!(condition_numbers, info.condition_number)
+            end
+            return cost_function(ψ, env′, operator)
+        end
+        g = only(gs)  # `withgradient` returns tuple of gradients `gs`
+        push!(gradnorms_unitcell, norm.(g.A))
+        push!(times, (time_ns() - start_time) * 1e-9)
+        return E, g
+    end
+
+    info = (
+        last_gradient=∂cost,
+        fg_evaluations=numfg,
+        costs=convergence_history[:, 1],
+        gradnorms=convergence_history[:, 2],
+        truncation_errors,
+        condition_numbers,
+        gradnorms_unitcell,
+        times,
+    )
+    return peps_final, env_final, cost, info
+end
+
+# Update PEPS unit cell in non-mutating way
+# Note: Both x and η are InfinitePEPS during optimization
+function peps_retract(x, η, α)
+    peps = deepcopy(x[1])
+    peps.A .+= η.A .* α
+    env = deepcopy(x[2])
+    return (peps, env), η
+end
+
+# Take real valued part of dot product
+real_inner(_, η₁, η₂) = real(dot(η₁, η₂))
+
+"""
+    symmetrize_retract_and_finalize!(symm::SymmetrizationStyle)
+
+Return the `retract` and `finalize!` function for symmetrizing the `peps` and `grad` tensors.
+"""
+function symmetrize_retract_and_finalize!(symm::SymmetrizationStyle)
+    finf = function symmetrize_finalize!((peps, env), E, grad, _)
+        grad_symm = symmetrize!(grad, symm)
+        return (peps, env), E, grad_symm
+    end
+    retf = function symmetrize_retract((peps, env), η, α)
+        peps_symm = deepcopy(peps)
+        peps_symm.A .+= η.A .* α
+        env′ = deepcopy(env)
+        symmetrize!(peps_symm, symm)
+        return (peps_symm, env′), η
+    end
+    return retf, finf
+end
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index 29a8f8287..7fdee6232 100644
--- a/src/algorithms/toolbox.jl
+++ b/src/algorithms/toolbox.jl
@@ -1,23 +1,23 @@
 """
-    expectation_value(peps::InfinitePEPS, O::LocalOperator, envs::CTMRGEnv)
+    expectation_value(peps::InfinitePEPS, O::LocalOperator, env::CTMRGEnv)
 
 Compute the expectation value ⟨peps|O|peps⟩ / ⟨peps|peps⟩ of a [`LocalOperator`](@ref) `O`
-for a PEPS `peps` using a given CTMRG environment `envs`.
+for a PEPS `peps` using a given CTMRG environment `env`.
 """
-function MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, envs::CTMRGEnv)
+function MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, env::CTMRGEnv)
     checklattice(peps, O)
     term_vals = dtmap([O.terms...]) do (inds, operator)  # OhMyThreads can't iterate over O.terms directly
-        contract_local_operator(inds, operator, peps, peps, envs) /
-        contract_local_norm(inds, peps, peps, envs)
+        contract_local_operator(inds, operator, peps, peps, env) /
+        contract_local_norm(inds, peps, peps, env)
     end
     return sum(term_vals)
 end
 """
-    expectation_value(pf::InfinitePartitionFunction, inds => O, envs::CTMRGEnv)
+    expectation_value(pf::InfinitePartitionFunction, inds => O, env::CTMRGEnv)
 
 Compute the expectation value corresponding to inserting a local tensor(s) `O` at
 position `inds` in the partition function `pf` and contracting the chole using a given CTMRG
-environment `envs`.
+environment `env`.
 
 Here `inds` can be specified as either a `Tuple{Int,Int}` or a `CartesianIndex{2}`, and `O`
 should be a rank-4 tensor conforming to the [`PartitionFunctionTensor`](@ref) indexing
@@ -26,19 +26,25 @@ convention.
 function MPSKit.expectation_value(
     pf::InfinitePartitionFunction,
     op::Pair{CartesianIndex{2},<:AbstractTensorMap{T,S,2,2}},
-    envs::CTMRGEnv,
+    env::CTMRGEnv,
 ) where {T,S}
-    return contract_local_tensor(op[1], op[2], envs) /
-           contract_local_tensor(op[1], pf[op[1]], envs)
+    return contract_local_tensor(op[1], op[2], env) /
+           contract_local_tensor(op[1], pf[op[1]], env)
 end
 function MPSKit.expectation_value(
-    pf::InfinitePartitionFunction, op::Pair{Tuple{Int,Int}}, envs::CTMRGEnv
+    pf::InfinitePartitionFunction, op::Pair{Tuple{Int,Int}}, env::CTMRGEnv
 )
-    return expectation_value(pf, CartesianIndex(op[1]) => op[2], envs)
+    return expectation_value(pf, CartesianIndex(op[1]) => op[2], env)
 end
 
-function costfun(peps::InfinitePEPS, envs::CTMRGEnv, O::LocalOperator)
-    E = MPSKit.expectation_value(peps, O, envs)
+"""
+    cost_function(peps::InfinitePEPS, env::CTMRGEnv, O::LocalOperator)
+
+Real part of expectation value of `O`. Prints a warning if the expectation value
+yields a finite imaginary part (up to a tolerance).
+"""
+function cost_function(peps::InfinitePEPS, env::CTMRGEnv, O::LocalOperator)
+    E = MPSKit.expectation_value(peps, O, env)
     ignore_derivatives() do
         isapprox(imag(E), 0; atol=sqrt(eps(real(E)))) ||
             @warn "Expectation value is not real: $E."
diff --git a/src/environments/vumps_environments.jl b/src/environments/vumps_environments.jl
index 1dd6ded89..b01f96f7f 100644
--- a/src/environments/vumps_environments.jl
+++ b/src/environments/vumps_environments.jl
@@ -4,18 +4,18 @@ using MPSKit: InfiniteEnvironments
 # ProductSpace instances
 
 function MPSKit.issamespace(
-    envs::InfiniteEnvironments,
+    env::InfiniteEnvironments,
     above::InfiniteMPS,
     operator::InfiniteTransferMatrix,
     below::InfiniteMPS,
 )
     L = MPSKit.check_length(above, operator, below)
     for i in 1:L
-        space(envs.GLs[i]) == (
+        space(env.GLs[i]) == (
             left_virtualspace(below, i) ⊗
             _elementwise_dual(left_virtualspace(operator, i)) ← left_virtualspace(above, i)
         ) || return false
-        space(envs.GRs[i]) == (
+        space(env.GRs[i]) == (
             right_virtualspace(above, i) ⊗ right_virtualspace(operator, i) ←
             right_virtualspace(below, i)
         ) || return false
diff --git a/src/operators/transfermatrix.jl b/src/operators/transfermatrix.jl
index abbf2496f..f39a19aea 100644
--- a/src/operators/transfermatrix.jl
+++ b/src/operators/transfermatrix.jl
@@ -231,10 +231,10 @@ the unit cell.
 
 @doc """
     MPSKit.leading_boundary(
-        st::InfiniteMPS, op::Union{InfiniteTransferPEPS,InfiniteTransferPEPO}, alg, [envs]
+        st::InfiniteMPS, op::Union{InfiniteTransferPEPS,InfiniteTransferPEPO}, alg, [env]
     )
     MPSKit.leading_boundary(
-        st::MPSMulitline, op::Union{MultilineTransferPEPS,MultilineTransferPEPO}, alg, [envs]
+        st::MPSMulitline, op::Union{MultilineTransferPEPS,MultilineTransferPEPO}, alg, [env]
     )
 
 Approximate the leading boundary MPS eigenvector for the transfer operator `op` using `st`
diff --git a/src/utility/symmetrization.jl b/src/utility/symmetrization.jl
index f1b519203..1f53f9303 100644
--- a/src/utility/symmetrization.jl
+++ b/src/utility/symmetrization.jl
@@ -216,23 +216,3 @@ function symmetrize!(peps::InfinitePEPS, symm::RotateReflect)
     end
     return peps
 end
-
-"""
-    symmetrize_retract_and_finalize!(symm::SymmetrizationStyle)
-
-Return the `retract` and `finalize!` function for symmetrizing the `peps` and `grad` tensors.
-"""
-function symmetrize_retract_and_finalize!(symm::SymmetrizationStyle)
-    finf = function symmetrize_finalize!((peps, envs), E, grad, _)
-        grad_symm = symmetrize!(grad, symm)
-        return (peps, envs), E, grad_symm
-    end
-    retf = function symmetrize_retract((peps, envs), η, α)
-        peps_symm = deepcopy(peps)
-        peps_symm.A .+= η.A .* α
-        e = deepcopy(envs)
-        symmetrize!(peps_symm, symm)
-        return (peps_symm, e), η
-    end
-    return retf, finf
-end
diff --git a/test/boundarymps/vumps.jl b/test/boundarymps/vumps.jl
index 272bc948e..8b816ad6a 100644
--- a/test/boundarymps/vumps.jl
+++ b/test/boundarymps/vumps.jl
@@ -14,10 +14,10 @@ const vumps_alg = VUMPS(; alg_eigsolve=MPSKit.Defaults.alg_eigsolve(; ishermitia
     T = PEPSKit.InfiniteTransferPEPS(psi, 1, 1)
     mps = PEPSKit.initializeMPS(T, [ComplexSpace(20)])
 
-    mps, envs, ϵ = leading_boundary(mps, T, vumps_alg)
+    mps, env, ϵ = leading_boundary(mps, T, vumps_alg)
     N = abs(sum(expectation_value(mps, T)))
 
-    ctm = leading_boundary(
+    ctm, = leading_boundary(
         CTMRGEnv(psi, ComplexSpace(20)), psi, SimultaneousCTMRG(; verbosity=1)
     )
     N´ = abs(norm(psi, ctm))
@@ -30,10 +30,10 @@ end
     T = PEPSKit.MultilineTransferPEPS(psi, 1)
 
     mps = PEPSKit.initializeMPS(T, fill(ComplexSpace(20), 2, 2))
-    mps, envs, ϵ = leading_boundary(mps, T, vumps_alg)
+    mps, env, ϵ = leading_boundary(mps, T, vumps_alg)
     N = abs(prod(expectation_value(mps, T)))
 
-    ctm = leading_boundary(
+    ctm, = leading_boundary(
         CTMRGEnv(psi, ComplexSpace(20)), psi, SimultaneousCTMRG(; verbosity=1)
     )
     N´ = abs(norm(psi, ctm))
@@ -68,6 +68,6 @@ end
     T = InfiniteTransferPEPO(psi, O, 1, 1)
 
     mps = PEPSKit.initializeMPS(T, [ComplexSpace(10)])
-    mps, envs, ϵ = leading_boundary(mps, T, vumps_alg)
+    mps, env, ϵ = leading_boundary(mps, T, vumps_alg)
     f = abs(prod(expectation_value(mps, T)))
 end
diff --git a/test/ctmrg/fixed_iterscheme.jl b/test/ctmrg/fixed_iterscheme.jl
index 208fdece0..0cf751837 100644
--- a/test/ctmrg/fixed_iterscheme.jl
+++ b/test/ctmrg/fixed_iterscheme.jl
@@ -30,7 +30,7 @@ atol = 1e-5
     # initialize states
     Random.seed!(2394823842)
     psi = InfinitePEPS(2, χbond; unitcell)
-    env_conv1 = leading_boundary(CTMRGEnv(psi, ComplexSpace(χenv)), psi, ctm_alg)
+    env_conv1, = leading_boundary(CTMRGEnv(psi, ComplexSpace(χenv)), psi, ctm_alg)
 
     # do extra iteration to get SVD
     env_conv2, info = ctmrg_iteration(psi, env_conv1, ctm_alg)
@@ -60,8 +60,8 @@ end
     # initialize states
     Random.seed!(91283219347)
     psi = InfinitePEPS(2, χbond)
-    env_init = CTMRGEnv(psi, ComplexSpace(χenv))
-    env_conv1 = leading_boundary(env_init, psi, ctm_alg_iter)
+    env₀ = CTMRGEnv(psi, ComplexSpace(χenv))
+    env_conv1, = leading_boundary(env₀, psi, ctm_alg_iter)
 
     # do extra iteration to get SVD
     env_conv2_iter, info_iter = ctmrg_iteration(psi, env_conv1, ctm_alg_iter)
diff --git a/test/ctmrg/flavors.jl b/test/ctmrg/flavors.jl
index 5d8b440c4..b3d935a54 100644
--- a/test/ctmrg/flavors.jl
+++ b/test/ctmrg/flavors.jl
@@ -19,10 +19,10 @@ projector_algs = [HalfInfiniteProjector, FullInfiniteProjector]
     # compute environments
     Random.seed!(32350283290358)
     psi = InfinitePEPS(2, χbond; unitcell)
-    env_sequential = leading_boundary(
+    env_sequential, = leading_boundary(
         CTMRGEnv(psi, ComplexSpace(χenv)), psi, ctm_alg_sequential
     )
-    env_simultaneous = leading_boundary(
+    env_simultaneous, = leading_boundary(
         CTMRGEnv(psi, ComplexSpace(χenv)), psi, ctm_alg_simultaneous
     )
 
@@ -50,8 +50,8 @@ projector_algs = [HalfInfiniteProjector, FullInfiniteProjector]
 
     # compare Heisenberg energies
     H = heisenberg_XYZ(InfiniteSquare(unitcell...))
-    E_sequential = costfun(psi, env_sequential, H)
-    E_simultaneous = costfun(psi, env_simultaneous, H)
+    E_sequential = cost_function(psi, env_sequential, H)
+    E_simultaneous = cost_function(psi, env_simultaneous, H)
     @test E_sequential ≈ E_simultaneous rtol = 1e-3
 end
 
@@ -68,7 +68,7 @@ end
     χs = [16 17 18; 15 20 21; 14 19 22]
     psi = InfinitePEPS(Ds, Ds, Ds)
     env = CTMRGEnv(psi, rand(10:20, 3, 3), rand(10:20, 3, 3))
-    env2 = leading_boundary(env, psi, ctm_alg)
+    env2, = leading_boundary(env, psi, ctm_alg)
 
     # check that the space is fixed
     @test all(space.(env.corners) .== space.(env2.corners))
diff --git a/test/ctmrg/gaugefix.jl b/test/ctmrg/gaugefix.jl
index 13d9a15b3..522314983 100644
--- a/test/ctmrg/gaugefix.jl
+++ b/test/ctmrg/gaugefix.jl
@@ -20,7 +20,7 @@ function _pre_converge_env(
     Random.seed!(seed)  # Seed RNG to make random environment consistent
     psi = InfinitePEPS(rand, T, physical_space, peps_space; unitcell)
     env₀ = CTMRGEnv(psi, ctm_space)
-    env_conv = leading_boundary(env₀, psi, SequentialCTMRG(; tol))
+    env_conv, = leading_boundary(env₀, psi, SequentialCTMRG(; tol))
     return env_conv, psi
 end
 
@@ -45,7 +45,7 @@ end
     alg = ctmrg_alg(; tol, projector_alg)
     env_pre, psi = preconv[(S, T, unitcell)]
     env_pre
-    env = leading_boundary(env_pre, psi, alg)
+    env, = leading_boundary(env_pre, psi, alg)
     env′, = ctmrg_iteration(psi, env, alg)
     env_fixed, = gauge_fix(env, env′)
     @test calc_elementwise_convergence(env, env_fixed) ≈ 0 atol = atol
diff --git a/test/ctmrg/gradients.jl b/test/ctmrg/gradients.jl
index 648d48bdb..6c55b5d63 100644
--- a/test/ctmrg/gradients.jl
+++ b/test/ctmrg/gradients.jl
@@ -50,28 +50,27 @@ steps = -0.01:0.005:0.01
     models
 )
     Pspace = Pspaces[i]
-    Vspace = Pspaces[i]
+    Vspace = Vspaces[i]
     Espace = Espaces[i]
     gms = gradmodes[i]
     calgs = ctmrg_algs[i]
-    psi_init = InfinitePEPS(Pspace, Vspace, Vspace)
     @testset "$ctmrg_alg and $alg_rrule" for (ctmrg_alg, alg_rrule) in
                                              Iterators.product(calgs, gms)
         @info "optimtest of $ctmrg_alg and $alg_rrule on $(names[i])"
         Random.seed!(42039482030)
-        dir = InfinitePEPS(Pspace, Vspace, Vspace)
-        psi = InfinitePEPS(Pspace, Vspace, Vspace)
-        env = leading_boundary(CTMRGEnv(psi, Espace), psi, ctmrg_alg)
+        dir = InfinitePEPS(Pspace, Vspace)
+        psi = InfinitePEPS(Pspace, Vspace)
+        env, = leading_boundary(CTMRGEnv(psi, Espace), psi, ctmrg_alg)
         alphas, fs, dfs1, dfs2 = OptimKit.optimtest(
             (psi, env),
             dir;
             alpha=steps,
             retract=PEPSKit.peps_retract,
             inner=PEPSKit.real_inner,
-        ) do (peps, envs)
+        ) do (peps, env)
             E, g = Zygote.withgradient(peps) do psi
-                envs2 = PEPSKit.hook_pullback(leading_boundary, envs, psi, ctmrg_alg; alg_rrule)
-                return costfun(psi, envs2, models[i])
+                env2, = PEPSKit.hook_pullback(leading_boundary, env, psi, ctmrg_alg; alg_rrule)
+                return cost_function(psi, env2, models[i])
             end
 
             return E, only(g)
diff --git a/test/ctmrg/jacobian_real_linear.jl b/test/ctmrg/jacobian_real_linear.jl
index d91c45294..24a655bfc 100644
--- a/test/ctmrg/jacobian_real_linear.jl
+++ b/test/ctmrg/jacobian_real_linear.jl
@@ -19,12 +19,12 @@ Dbond, χenv = 2, 16
 @testset "$iterscheme and $ctm_alg" for (iterscheme, ctm_alg) in algs
     Random.seed!(123521938519)
     state = InfinitePEPS(2, Dbond)
-    envs = leading_boundary(CTMRGEnv(state, ComplexSpace(χenv)), state, ctm_alg)
+    env, = leading_boundary(CTMRGEnv(state, ComplexSpace(χenv)), state, ctm_alg)
 
     # follow code of _rrule
     if iterscheme == :fixed
-        envs_conv, info = ctmrg_iteration(state, envs, ctm_alg)
-        envs_fixed, signs = gauge_fix(envs, envs_conv)
+        env_conv, info = ctmrg_iteration(state, env, ctm_alg)
+        env_fixed, signs = gauge_fix(env, env_conv)
         U_fixed, V_fixed = fix_relative_phases(info.U, info.V, signs)
         svd_alg_fixed = SVDAdjoint(;
             fwd_alg=PEPSKit.FixedSVD(U_fixed, info.S, V_fixed),
@@ -33,19 +33,19 @@ Dbond, χenv = 2, 16
         alg_fixed = @set ctm_alg.projector_alg.svd_alg = svd_alg_fixed
         alg_fixed = @set alg_fixed.projector_alg.trscheme = notrunc()
 
-        _, env_vjp = pullback(state, envs_fixed) do A, x
+        _, env_vjp = pullback(state, env_fixed) do A, x
             e, = PEPSKit.ctmrg_iteration(A, x, alg_fixed)
             return PEPSKit.fix_global_phases(x, e)
         end
     elseif iterscheme == :diffgauge
-        _, env_vjp = pullback(state, envs) do A, x
+        _, env_vjp = pullback(state, env) do A, x
             return gauge_fix(x, ctmrg_iteration(A, x, ctm_alg)[1])[1]
         end
     end
 
     # get Jacobians of single iteration
     ∂f∂A(x)::typeof(state) = env_vjp(x)[1]
-    ∂f∂x(x)::typeof(envs) = env_vjp(x)[2]
+    ∂f∂x(x)::typeof(env) = env_vjp(x)[2]
 
     # compute real and complex errors
     env_in = CTMRGEnv(state, ComplexSpace(16))
diff --git a/test/ctmrg/partition_function.jl b/test/ctmrg/partition_function.jl
index 8995d6ab8..13a20ffb1 100644
--- a/test/ctmrg/partition_function.jl
+++ b/test/ctmrg/partition_function.jl
@@ -98,7 +98,7 @@ projector_algs = [HalfInfiniteProjector, FullInfiniteProjector]
     ctm_styles, projector_algs
 )
     ctm_alg = ctm_style(; maxiter=150, projector_alg)
-    env = leading_boundary(env0, Z, ctm_alg)
+    env, = leading_boundary(env0, Z, ctm_alg)
 
     # check observables
     λ = PEPSKit.value(Z, env)
diff --git a/test/heisenberg.jl b/test/heisenberg.jl
index 48e58be01..24ea65ae9 100644
--- a/test/heisenberg.jl
+++ b/test/heisenberg.jl
@@ -21,14 +21,14 @@ E_ref = -0.6602310934799577
     # initialize states
     Random.seed!(123)
     H = heisenberg_XYZ(InfiniteSquare())
-    psi_init = InfinitePEPS(2, Dbond)
-    env_init = leading_boundary(CTMRGEnv(psi_init, ComplexSpace(χenv)), psi_init, ctm_alg)
+    peps₀ = InfinitePEPS(2, Dbond)
+    env₀, = leading_boundary(CTMRGEnv(peps₀, ComplexSpace(χenv)), peps₀, ctm_alg)
 
     # optimize energy and compute correlation lengths
-    result = fixedpoint(psi_init, H, opt_alg, env_init)
-    ξ_h, ξ_v, = correlation_length(result.peps, result.env)
+    peps, env, E, = fixedpoint(H, peps₀, env₀, opt_alg)
+    ξ_h, ξ_v, = correlation_length(peps, env)
 
-    @test result.E ≈ E_ref atol = 1e-2
+    @test E ≈ E_ref atol = 1e-2
     @test all(@. ξ_h > 0 && ξ_v > 0)
 end
 
@@ -36,18 +36,16 @@ end
     # initialize states
     Random.seed!(456)
     unitcell = (1, 2)
-    H_1x2 = heisenberg_XYZ(InfiniteSquare(unitcell...))
-    psi_init_1x2 = InfinitePEPS(2, Dbond; unitcell)
-    env_init_1x2 = leading_boundary(
-        CTMRGEnv(psi_init_1x2, ComplexSpace(χenv)), psi_init_1x2, ctm_alg
-    )
+    H = heisenberg_XYZ(InfiniteSquare(unitcell...))
+    peps₀ = InfinitePEPS(2, Dbond; unitcell)
+    env₀, = leading_boundary(CTMRGEnv(peps₀, ComplexSpace(χenv)), peps₀, ctm_alg)
 
     # optimize energy and compute correlation lengths
-    result_1x2 = fixedpoint(psi_init_1x2, H_1x2, opt_alg, env_init_1x2)
-    ξ_h_1x2, ξ_v_1x2, = correlation_length(result_1x2.peps, result_1x2.env)
+    peps, env, E, = fixedpoint(H, peps₀, env₀, opt_alg)
+    ξ_h, ξ_v, = correlation_length(peps, env)
 
-    @test result_1x2.E ≈ 2 * E_ref atol = 1e-2
-    @test all(@. ξ_h_1x2 > 0 && ξ_v_1x2 > 0)
+    @test E ≈ 2 * E_ref atol = 1e-2
+    @test all(@. ξ_h > 0 && ξ_v > 0)
 end
 
 @testset "Simple update into AD optimization" begin
@@ -82,11 +80,11 @@ end
 
     # absorb weight into site tensors and CTMRG
     peps = InfinitePEPS(wpeps)
-    envs₀ = CTMRGEnv(rand, Float64, peps, Espace)
-    envs = leading_boundary(envs₀, peps, SimultaneousCTMRG())
+    env₀ = CTMRGEnv(rand, Float64, peps, Espace)
+    env, = leading_boundary(env₀, peps, SimultaneousCTMRG())
 
     # measure physical quantities
-    e_site = costfun(peps, envs, ham) / (N1 * N2)
+    e_site = cost_function(peps, env, ham) / (N1 * N2)
     @info "Simple update energy = $e_site"
     # benchmark data from Phys. Rev. B 94, 035133 (2016)
     @test isapprox(e_site, -0.6594; atol=1e-3)
@@ -94,9 +92,9 @@ end
     # continue with auto differentiation
     svd_alg_gmres = SVDAdjoint(; rrule_alg=GMRES(; tol=1e-5))
     opt_alg_gmres = @set opt_alg.boundary_alg.projector_alg.svd_alg = svd_alg_gmres
-    result = fixedpoint(peps, ham, opt_alg_gmres, envs)  # sensitivity warnings and degeneracies due to SU(2)?
-    ξ_h, ξ_v, = correlation_length(result.peps, result.env)
-    e_site2 = result.E / (N1 * N2)
+    peps_final, env_final, E_final, = fixedpoint(ham, peps, env, opt_alg_gmres)  # sensitivity warnings and degeneracies due to SU(2)?
+    ξ_h, ξ_v, = correlation_length(peps_final, env_final)
+    e_site2 = E_final / (N1 * N2)
     @info "Auto diff energy = $e_site2"
     @test e_site2 ≈ E_ref atol = 1e-2
     @test all(@. ξ_h > 0 && ξ_v > 0)
diff --git a/test/j1j2_model.jl b/test/j1j2_model.jl
index 5d13bba74..391f81a21 100644
--- a/test/j1j2_model.jl
+++ b/test/j1j2_model.jl
@@ -13,22 +13,23 @@ opt_alg = PEPSOptimize(;
     boundary_alg=ctm_alg,
     optimizer=LBFGS(4; gradtol=1e-3, verbosity=3),
     gradient_alg=LinSolver(; iterscheme=:diffgauge),
+    symmetrization=RotateReflect(),
 )
 
 # initialize states
 Random.seed!(91283219347)
 H = j1_j2(InfiniteSquare(); J2=0.25)
-psi_init = product_peps(2, χbond; noise_amp=1e-1)
-psi_init = symmetrize!(psi_init, RotateReflect())
-env_init = leading_boundary(CTMRGEnv(psi_init, ComplexSpace(χenv)), psi_init, ctm_alg);
+peps₀ = product_peps(2, χbond; noise_amp=1e-1)
+peps₀ = symmetrize!(peps₀, RotateReflect())
+env₀, = leading_boundary(CTMRGEnv(peps₀, ComplexSpace(χenv)), peps₀, ctm_alg);
 
 # find fixedpoint
-result = fixedpoint(psi_init, H, opt_alg, env_init; symmetrization=RotateReflect())
-ξ_h, ξ_v, = correlation_length(result.peps, result.env)
+peps, env, E, = fixedpoint(H, peps₀, env₀, opt_alg)
+ξ_h, ξ_v, = correlation_length(peps, env)
 
 # compare against Juraj Hasik's data:
 # https://github.com/jurajHasik/j1j2_ipeps_states/blob/main/single-site_pg-C4v-A1/j20.25/state_1s_A1_j20.25_D2_chi_opt48.dat
 ξ_ref = -1 / log(0.2723596743547324)
-@test result.E ≈ -0.5618837021945925 atol = 1e-3
+@test E ≈ -0.5618837021945925 atol = 1e-3
 @test all(@. isapprox(ξ_h, ξ_ref; atol=1e-1) && isapprox(ξ_v, ξ_ref; atol=1e-1))
 @test ξ_h ≈ ξ_v atol = 1e-6  # Test symmetrization of optimized PEPS and environment
diff --git a/test/pwave.jl b/test/pwave.jl
index 4e36b14e9..00a65391a 100644
--- a/test/pwave.jl
+++ b/test/pwave.jl
@@ -8,23 +8,23 @@ using OptimKit
 # Initialize parameters
 unitcell = (2, 2)
 H = pwave_superconductor(InfiniteSquare(unitcell...))
-χbond = 2
+Dbond = 2
 χenv = 16
-ctm_alg = SimultaneousCTMRG(; maxiter=150)
+ctm_alg = SimultaneousCTMRG()
 opt_alg = PEPSOptimize(;
     boundary_alg=ctm_alg, optimizer=LBFGS(4; maxiter=10, gradtol=1e-3, verbosity=3)
 )
 
 # initialize states
 Pspace = Vect[FermionParity](0 => 1, 1 => 1)
-Vspace = Vect[FermionParity](0 => χbond ÷ 2, 1 => χbond ÷ 2)
+Vspace = Vect[FermionParity](0 => Dbond ÷ 2, 1 => Dbond ÷ 2)
 Envspace = Vect[FermionParity](0 => χenv ÷ 2, 1 => χenv ÷ 2)
 Random.seed!(91283219347)
-psi_init = InfinitePEPS(Pspace, Vspace, Vspace; unitcell)
-env_init = leading_boundary(CTMRGEnv(psi_init, Envspace), psi_init, ctm_alg);
+peps₀ = InfinitePEPS(Pspace, Vspace, Vspace; unitcell)
+env₀, = leading_boundary(CTMRGEnv(peps₀, Envspace), peps₀, ctm_alg);
 
 # find fixedpoint
-result = fixedpoint(psi_init, H, opt_alg, env_init)
+_, _, E, = fixedpoint(H, peps₀, env₀, opt_alg)
 
 # comparison with Gaussian PEPS minimum at D=2 on 1000x1000 square lattice with aPBC
-@test result.E / prod(size(psi_init)) ≈ -2.6241 atol = 5e-2
+@test E / prod(size(peps₀)) ≈ -2.6241 atol = 5e-2
diff --git a/test/tf_ising.jl b/test/tf_ising.jl
index 903770214..f4ad251c0 100644
--- a/test/tf_ising.jl
+++ b/test/tf_ising.jl
@@ -27,25 +27,25 @@ opt_alg = PEPSOptimize(;
 # initialize states
 H = transverse_field_ising(InfiniteSquare(); g)
 Random.seed!(2928528935)
-psi_init = InfinitePEPS(2, χbond)
-env_init = leading_boundary(CTMRGEnv(psi_init, ComplexSpace(χenv)), psi_init, ctm_alg)
+peps₀ = InfinitePEPS(2, χbond)
+env₀, = leading_boundary(CTMRGEnv(peps₀, ComplexSpace(χenv)), peps₀, ctm_alg)
 
 # find fixedpoint
-result = fixedpoint(psi_init, H, opt_alg, env_init)
-ξ_h, ξ_v, = correlation_length(result.peps, result.env)
+peps, env, E, = fixedpoint(H, peps₀, env₀, opt_alg)
+ξ_h, ξ_v, = correlation_length(peps, env)
 
 # compute magnetization
-σx = TensorMap(scalartype(psi_init)[0 1; 1 0], ℂ^2, ℂ^2)
+σx = TensorMap(scalartype(peps₀)[0 1; 1 0], ℂ^2, ℂ^2)
 M = LocalOperator(H.lattice, (CartesianIndex(1, 1),) => σx)
-magn = expectation_value(result.peps, M, result.env)
+magn = expectation_value(peps, M, env)
 
-@test result.E ≈ e atol = 1e-2
+@test E ≈ e atol = 1e-2
 @test imag(magn) ≈ 0 atol = 1e-6
 @test abs(magn) ≈ mˣ atol = 5e-2
 
 # find fixedpoint in polarized phase and compute correlations lengths
 H_polar = transverse_field_ising(InfiniteSquare(); g=4.5)
-result_polar = fixedpoint(psi_init, H_polar, opt_alg, env_init)
-ξ_h_polar, ξ_v_polar, = correlation_length(result_polar.peps, result_polar.env)
+peps_polar, env_polar, = fixedpoint(H_polar, peps₀, env₀, opt_alg)
+ξ_h_polar, ξ_v_polar, = correlation_length(peps_polar, env_polar)
 @test ξ_h_polar < ξ_h
 @test ξ_v_polar < ξ_v