diff --git a/Project.toml b/Project.toml
index 0dfbe5e4d..c73ef4802 100644
--- a/Project.toml
+++ b/Project.toml
@@ -45,7 +45,7 @@ OptimKit = "0.5"
 Printf = "1"
 Random = "1"
 Statistics = "1"
-TensorKit = "0.16.5, 0.17"
+TensorKit = "0.17"
 TensorKitTensors = "0.3.1"
 TensorOperations = "5"
 TupleTools = "1.6.0"
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 75e876381..7196cfa2a 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -32,6 +32,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 
 ### Fixed
 
+- `initializePEPS` now returns a real-valued PEPS if given a real-valued input PEPO, rather than always returning a complex-valued PEPS
+
 ### Performance
 
 ## [0.8.1](https://github.com/quantumkithub/pepskit.jl/compare/v0.8.0...v0.8.1) - 2026-07-30
diff --git a/src/algorithms/ctmrg/gaugefix.jl b/src/algorithms/ctmrg/gaugefix.jl
index 6ff44e554..9c5c73d78 100644
--- a/src/algorithms/ctmrg/gaugefix.jl
+++ b/src/algorithms/ctmrg/gaugefix.jl
@@ -82,7 +82,7 @@ function compute_relative_phases(
 
         # Random MPS of same bond dimension
         M = map(Tsfinal) do t
-            randn(scalartype(t), codomain(t) ← domain(t))
+            randn(storagetype(T), codomain(t) ← domain(t))
         end
 
         # Find right fixed points of mixed transfer matrices
@@ -107,7 +107,7 @@ function compute_relative_phases(envfinal::CTMRGEnv{C, T}, envprev::CTMRGEnv{C,
 
     # Random Hermitian MPS of same bond dimension
     # (make Hermitian such that T-M transfer matrix has real eigenvalues)
-    M = _project_hermitian(randn(scalartype(Tfinal), space(Tfinal)))
+    M = _project_hermitian(randn(storagetype(Tfinal), space(Tfinal)))
 
     # Find right fixed points of mixed transfer matrices
     eigsolve_alg = Lanczos() # real eigenvalues
@@ -145,7 +145,7 @@ end
 
 function initialize_right_fixedpoint(tops, bottoms)
     ρ0 = randn(
-        scalartype(tops), space(tops[end], numind(tops[end]))' ← space(bottoms[end], numind(bottoms[end]))'
+        TensorKit.promote_storagetype(tops...), space(tops[end], numind(tops[end]))' ← space(bottoms[end], numind(bottoms[end]))'
     )
     return ρ0
 end
diff --git a/src/algorithms/ctmrg/initialization.jl b/src/algorithms/ctmrg/initialization.jl
index 235611e0b..42d7a6d87 100644
--- a/src/algorithms/ctmrg/initialization.jl
+++ b/src/algorithms/ctmrg/initialization.jl
@@ -1,11 +1,11 @@
 """
-    initialize_ctmrg_environment([elt::Type{<:Number},] n::InfiniteSquareNetwork, alg::RandomInitialization, virtual_spaces...)
+    initialize_ctmrg_environment([elt::Type,] n::InfiniteSquareNetwork, alg::RandomInitialization, virtual_spaces...)
 
 Initialize a fully random `CTMRGEnv` using the given environment virtual spaces. See
 [`CTMRGEnv`](@ref) for details on the expected format of the virtual spaces.
 """
 function initialize_ctmrg_environment(
-        elt::Type{<:Number},
+        elt::Type,
         n::InfiniteSquareNetwork,
         alg::RandomInitialization,
         virtual_spaces... = oneunit(spacetype(n)),
@@ -14,13 +14,13 @@ function initialize_ctmrg_environment(
 end
 
 """
-    initialize_ctmrg_environment([elt::Type{<:Number},] n::InfiniteSquareNetwork, alg::RandomInitialization)
+    initialize_ctmrg_environment([elt::Type,] n::InfiniteSquareNetwork, alg::RandomInitialization)
 
 Initialize a `CTMRGEnv` corresponding to a product state with trivial virtual spaces and
 corners. The product state edge tensors are initialized as `alg.f(elt, V::ProductSpace)`.
 """
 function initialize_ctmrg_environment(
-        elt::Type{<:Number},
+        elt::Type,
         n::InfiniteSquareNetwork,
         alg::ProductStateInitialization,
     )
@@ -29,14 +29,14 @@ function initialize_ctmrg_environment(
 end
 
 """
-    initialize_ctmrg_environment([elt::Type{<:Number},] n::InfiniteSquareNetwork, alg::RandomInitialization, [env0])
+    initialize_ctmrg_environment([elt::Type,] n::InfiniteSquareNetwork, alg::RandomInitialization, [env0])
 
 Initialize a `CTMRGEnv` by applying a single untruncated iteration of
 [`SimultaneousCTMRG`](@ref) to a given initial environment. By default, the starting
 environment is chosen as a random product state.
 """
 function initialize_ctmrg_environment(
-        elt::Type{<:Number},
+        elt::Type,
         n::InfiniteSquareNetwork,
         alg::ApplicationInitialization,
         env0 = ProductStateEnv(alg.f, elt, n)
@@ -50,7 +50,7 @@ _check_two_layer(::InfiniteSquareNetwork) = false
 _check_two_layer(::InfiniteSquareNetwork{<:PEPSSandwich}) = true
 
 """
-    initialize_ctmrg_environment([elt::Type{<:Number},] n::InfiniteSquareNetwork, alg::RandomInitialization, [env0])
+    initialize_ctmrg_environment([elt::Type,] n::InfiniteSquareNetwork, alg::RandomInitialization, [env0])
 
 Initialize a `CTMRGEnv` corresponding to a product state acting as an identity between the
 virtual spaces of a two-layer network, for example
@@ -65,7 +65,7 @@ virtual spaces of a two-layer network, for example
 ```
 """
 function initialize_ctmrg_environment(
-        elt::Type{<:Number},
+        elt::Type,
         n::InfiniteSquareNetwork,
         ::IdentityInitialization,
     )
@@ -86,5 +86,5 @@ function initialize_ctmrg_environment(
         elt::Type{<:Number}, A::Union{InfinitePEPS, InfinitePartitionFunction}, args...;
         kwargs...
     )
-    return initialize_ctmrg_environment(elt, InfiniteSquareNetwork(A), args...; kwargs...)
+    return initialize_ctmrg_environment(similarstoragetype(storagetype(A), elt), InfiniteSquareNetwork(A), args...; kwargs...)
 end
diff --git a/src/algorithms/time_evolution/gaugefix_su.jl b/src/algorithms/time_evolution/gaugefix_su.jl
index 13a2782d6..a373d773d 100644
--- a/src/algorithms/time_evolution/gaugefix_su.jl
+++ b/src/algorithms/time_evolution/gaugefix_su.jl
@@ -17,7 +17,7 @@ $(TYPEDFIELDS)
     maxiter::Int = 100
 end
 
-function _trivial_gates(elt::Type{<:Number}, lattice::Matrix{S}) where {S <: ElementarySpace}
+function _trivial_gates(elt::Type, lattice::Matrix{S}) where {S <: ElementarySpace}
     Nr, Nc = size(lattice)
     gates = map(Iterators.product(1:2, 1:Nc, 1:Nr)) do (d, c, r)
         site1 = CartesianIndex(r, c)
@@ -38,7 +38,7 @@ Fix the gauge of `psi` using trivial simple update.
 """
 function gauge_fix(psi::InfiniteState, alg::SUGauge)
     time0 = time()
-    gates = _trivial_gates(scalartype(psi), physicalspace(psi))
+    gates = _trivial_gates(storagetype(psi), physicalspace(psi))
     trunc = _get_fixedspacetrunc(psi)
     su_alg = SimpleUpdate(; trunc, bipartite = _is_bipartite(psi))
     wts0 = SUWeight(psi)
diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl
index da79055e2..2b5b87673 100644
--- a/src/algorithms/time_evolution/simpleupdate3site.jl
+++ b/src/algorithms/time_evolution/simpleupdate3site.jl
@@ -1,6 +1,6 @@
 function _fuse_physicalspaces(O::GenericMPSTensor{S, 5}) where {S <: ElementarySpace}
     V1, V2 = codomain(O, 2), codomain(O, 3)
-    F = isomorphism(Int, fuse(V1, V2), V1 ⊗ V2)
+    F = isomorphism(similarstoragetype(O, Int), fuse(V1, V2), V1 ⊗ V2)
     @plansor O_fused[-1 -2 -4 -5; -6] := F[-2; 2 3] * O[-1 2 3 -4 -5; -6]
     return O_fused, F
 end
@@ -8,7 +8,7 @@ end
 function _unfuse_physicalspace(
         O::GenericMPSTensor{S, 4}, Vout::ElementarySpace, Vin::ElementarySpace = Vout'
     ) where {S <: ElementarySpace}
-    F = isomorphism(Int, Vout ⊗ Vin, fuse(Vout ⊗ Vin))
+    F = isomorphism(similarstoragetype(O, Int), Vout ⊗ Vin, fuse(Vout ⊗ Vin))
     @plansor O_unfused[-1 -2 -3 -4 -5; -6] := F[-2 -3; 1] * O[-1 1 -4 -5; -6]
     return O_unfused, F
 end
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
index c1dc32233..d6065d20a 100644
--- a/src/algorithms/time_evolution/time_evolve.jl
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -69,7 +69,7 @@ function _timeevol_sanity_check(
 end
 
 function MPSKit.infinite_temperature_density_matrix(H::LocalOperator)
-    T = scalartype(H)
+    T = storagetype(H)
     A = map(physicalspace(H)) do Vp
         ψ = permute(TensorKit.id(T, Vp), (1, 2))
         Vv = oneunit(Vp) # trivial (1D) virtual space
diff --git a/src/algorithms/time_evolution/trotter_gate.jl b/src/algorithms/time_evolution/trotter_gate.jl
index c4a52bfd5..9503ca4d2 100644
--- a/src/algorithms/time_evolution/trotter_gate.jl
+++ b/src/algorithms/time_evolution/trotter_gate.jl
@@ -157,7 +157,9 @@ function _trotterize_nnn2site!(gates::Vector, H::LocalOperator, dt::Number)
             term = permute(term, ((2, 1), (4, 3)))
         end
         gate = gate_to_mpo(exp(term * -dt / 2))
-        b = TensorKit.BraidingTensor{T}(physicalspace(H, x2), left_virtualspace(gate[2]))
+        A = similarstoragetype(term, T)
+        S = spacetype(H)
+        b = TensorKit.BraidingTensor{T, S, A}(physicalspace(H, x2), left_virtualspace(gate[2]))
         insert!(gate, 2, TensorMap(b))
         push!(gates, [x1, x2, x3] => gate)
     end
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index 533fac679..b825b2f80 100644
--- a/src/algorithms/toolbox.jl
+++ b/src/algorithms/toolbox.jl
@@ -13,7 +13,7 @@ function edge_transfer_spectrum(
         sector = one(sectortype(E))
     ) where {E <: CTMRGEdgeTensor}
     init = randn(
-        scalartype(E),
+        storagetype(E),
         space(first(bot), numind(first(bot)))' ← ℂ[typeof(sector)](sector => 1)' ⊗ space(first(top), 1),
     )
 
diff --git a/src/environments/bp_environments.jl b/src/environments/bp_environments.jl
index 70460f721..afc2856f9 100644
--- a/src/environments/bp_environments.jl
+++ b/src/environments/bp_environments.jl
@@ -108,7 +108,7 @@ function BPEnv(
         f, T, D_north::P, D_east::P;
         unitcell::Tuple{Int, Int} = (1, 1), posdef::Bool = true
     ) where {P <: ProductSpace}
-    return BPEnv(f, T, N, _fill_edge_physical_spaces(D_north, D_east; unitcell)...; posdef)
+    return BPEnv(f, T, _fill_edge_physical_spaces(D_north, D_east; unitcell)...; posdef)
 end
 
 """
@@ -119,13 +119,13 @@ Construct a BP environment by specifying a corresponding [`InfiniteSquareNetwork
 function BPEnv(f, T, network::InfiniteSquareNetwork; posdef::Bool = true)
     Ds_north = _north_edge_physical_spaces(network)
     Ds_east = _east_edge_physical_spaces(network)
-    return BPEnv(f, T, Ds_north, Ds_east; posdef)
+    return BPEnv(f, similarstoragetype(storagetype(network), eltype(T)), Ds_north, Ds_east; posdef)
 end
 function BPEnv(network::Union{InfiniteSquareNetwork, InfinitePartitionFunction, InfinitePEPS, InfinitePEPO}, args...; kwargs...)
-    return BPEnv(isomorphism, scalartype(network), network, args...; kwargs...)
+    return BPEnv(isomorphism, storagetype(network), network, args...; kwargs...)
 end
 function BPEnv(f, T, state::Union{InfinitePartitionFunction, InfinitePEPS, InfinitePEPO}, args...; kwargs...)
-    return BPEnv(f, T, InfiniteSquareNetwork(state), args...; kwargs...)
+    return BPEnv(f, similarstoragetype(eltype(state), eltype(T)), InfiniteSquareNetwork(state), args...; kwargs...)
 end
 
 Base.eltype(::Type{BPEnv{T}}) where {T} = T
@@ -173,7 +173,7 @@ function CTMRGEnv(bp_env::BPEnv)
         return insertleftunit(insertleftunit(M), 1)
     end
     corners = map(CartesianIndices(edges)) do _
-        return TensorKit.id(scalartype(bp_env), oneunit(spacetype(bp_env)))
+        return TensorKit.id(storagetype(bp_env), oneunit(spacetype(bp_env)))
     end
     return CTMRGEnv(corners, edges)
 end
diff --git a/src/environments/product_state_environments.jl b/src/environments/product_state_environments.jl
index 394bc5751..93a1577b9 100644
--- a/src/environments/product_state_environments.jl
+++ b/src/environments/product_state_environments.jl
@@ -82,13 +82,13 @@ Construct a product state environment by specifying a corresponding [`InfiniteSq
 function ProductStateEnv(f, T, network::InfiniteSquareNetwork)
     Ds_north = _north_edge_physical_spaces(network)
     Ds_east = _east_edge_physical_spaces(network)
-    return ProductStateEnv(f, T, Ds_north, Ds_east)
+    return ProductStateEnv(f, similarstoragetype(storagetype(network), eltype(T)), Ds_north, Ds_east)
 end
 function ProductStateEnv(network::Union{InfiniteSquareNetwork, InfinitePartitionFunction, InfinitePEPS})
-    return ProductStateEnv(randn, scalartype(network), network)
+    return ProductStateEnv(randn, storagetype(network), network)
 end
 function ProductStateEnv(f, T, state::Union{InfinitePartitionFunction, InfinitePEPS}, args...)
-    return ProductStateEnv(f, T, InfiniteSquareNetwork(state), args...)
+    return ProductStateEnv(f, similarstoragetype(eltype(state), eltype(T)), InfiniteSquareNetwork(state), args...)
 end
 
 Base.eltype(::Type{ProductStateEnv{T}}) where {T} = T
diff --git a/src/environments/suweight.jl b/src/environments/suweight.jl
index 3ccd5ba38..796c4e711 100644
--- a/src/environments/suweight.jl
+++ b/src/environments/suweight.jl
@@ -58,21 +58,43 @@ function SUWeight(wts_mats::AbstractMatrix{E}...) where {E <: PEPSWeight}
 end
 
 """
-    SUWeight(Nspaces::M, [Espaces::M]) where {M<:AbstractMatrix{<:ElementarySpace}}
+    SUWeight([T = ComplexF64,] Nspaces::M, [Espaces::M]) where {M<:AbstractMatrix{<:ElementarySpace}}
 
 Create a trivial `SUWeight` by specifying the vertical (north) or horizontal (east) virtual bond spaces.
+`T` is either a scalar type (like `ComplexF32`) or an array type (like `Vector{Float64}`) that
+specifies the storage type. If `T <: Number` the default storage type is `Vector{T}`.
 """
 function SUWeight(
+        T,
         Nspaces::M, Espaces::M = Nspaces
     ) where {M <: AbstractMatrix{<:ElementarySpace}}
     @assert size(Nspaces) == size(Espaces)
     Nr, Nc = size(Nspaces)
     weights = map(Iterators.product(1:2, 1:Nr, 1:Nc)) do (d, r, c)
         V = (d == 1 ? Espaces[r, c] : Nspaces[r, c])
-        DiagonalTensorMap(ones(reduceddim(V)), V)
+        realT = similarstoragetype(similarstoragetype(T), real(eltype(T)))
+        diag = realT(undef, reduceddim(V))
+        fill!(diag, one(real(eltype(T))))
+        DiagonalTensorMap(diag, V)
     end
     return SUWeight(weights)
 end
+SUWeight(Nspaces::M, Espaces::M = Nspaces) where {M <: AbstractMatrix{<:ElementarySpace}} =
+    SUWeight(ComplexF64, Nspaces, Espaces)
+
+"""
+    SUWeight(T, Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int}=(1, 1)) where {S<:ElementarySpace}
+
+Create a trivial `SUWeight` by specifying its vertical (north) and horizontal (east)
+as `ElementarySpace`s) and unit cell size. `T` is either a scalar type (like `ComplexF32`)
+or an array type (like `Vector{Float64}`) that specifies the storage type.
+If `T <: Number` the default storage type is `Vector{T}`.
+"""
+function SUWeight(
+        T, Nspace::S, Espace::S = Nspace; unitcell::Tuple{Int, Int} = (1, 1)
+    ) where {S <: ElementarySpace}
+    return SUWeight(T, fill(Nspace, unitcell), fill(Espace, unitcell))
+end
 
 """
     SUWeight(Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int}=(1, 1)) where {S<:ElementarySpace}
@@ -83,7 +105,7 @@ as `ElementarySpace`s) and unit cell size.
 function SUWeight(
         Nspace::S, Espace::S = Nspace; unitcell::Tuple{Int, Int} = (1, 1)
     ) where {S <: ElementarySpace}
-    return SUWeight(fill(Nspace, unitcell), fill(Espace, unitcell))
+    return SUWeight(ComplexF64, fill(Nspace, unitcell), fill(Espace, unitcell))
 end
 
 """
@@ -94,7 +116,7 @@ Create a trivial `SUWeight` for a given InfinitePEPS.
 function SUWeight(peps::InfinitePEPS)
     Nspaces = map(Base.Fix2(domain, NORTH), unitcell(peps))
     Espaces = map(Base.Fix2(domain, EAST), unitcell(peps))
-    return SUWeight(Nspaces, Espaces)
+    return SUWeight(storagetype(peps), Nspaces, Espaces)
 end
 
 """
@@ -106,7 +128,7 @@ function SUWeight(pepo::InfinitePEPO)
     @assert size(pepo, 3) == 1
     Nspaces = map(Base.Fix2(domain, NORTH), @view(unitcell(pepo)[:, :, 1]))
     Espaces = map(Base.Fix2(domain, EAST), @view(unitcell(pepo)[:, :, 1]))
-    return SUWeight(Nspaces, Espaces)
+    return SUWeight(storagetype(pepo), Nspaces, Espaces)
 end
 
 Random.rand!(wts::SUWeight) = rand!(Random.default_rng(), wts)
@@ -329,7 +351,7 @@ which has the same real scalartype as ``wts`.
 """
 function CTMRGEnv(wts::SUWeight)
     _, Nr, Nc = size(wts)
-    elt = scalartype(wts)
+    elt = storagetype(wts)
     V_env = oneunit(spacetype(wts))
     edges = map(Iterators.product(1:4, 1:Nr, 1:Nc)) do (d, r, c)
         wt_idx = if d == NORTH
diff --git a/src/environments/vumps_environments.jl b/src/environments/vumps_environments.jl
index 48605f2f1..6ef0375d0 100644
--- a/src/environments/vumps_environments.jl
+++ b/src/environments/vumps_environments.jl
@@ -25,10 +25,11 @@ function MPSKit.allocate_GL(
         bra::InfiniteMPS, mpo::InfiniteTransferMatrix, ket::InfiniteMPS, i::Int
     )
     T = Base.promote_type(scalartype(bra), scalartype(mpo), scalartype(ket))
+    TA = similarstoragetype(storagetype(mpo), T)
     V =
         left_virtualspace(bra, i) ⊗ _elementwise_dual(left_virtualspace(mpo, i)) ←
         left_virtualspace(ket, i)
-    TT = TensorMap{T}
+    TT = TensorKit.TensorMapWithStorage{T, TA}
     return TT(undef, V)
 end
 
@@ -36,7 +37,8 @@ function MPSKit.allocate_GR(
         bra::InfiniteMPS, mpo::InfiniteTransferMatrix, ket::InfiniteMPS, i::Int
     )
     T = Base.promote_type(scalartype(bra), scalartype(mpo), scalartype(ket))
+    TA = similarstoragetype(storagetype(mpo), T)
     V = right_virtualspace(ket, i) ⊗ right_virtualspace(mpo, i) ← right_virtualspace(bra, i)
-    TT = TensorMap{T}
+    TT = TensorKit.TensorMapWithStorage{T, TA}
     return TT(undef, V)
 end
diff --git a/src/networks/tensors.jl b/src/networks/tensors.jl
index 6ac6411da..347860366 100644
--- a/src/networks/tensors.jl
+++ b/src/networks/tensors.jl
@@ -24,10 +24,11 @@ const PFTensor = PartitionFunctionTensor
 
 """
     PartitionFunctionTensor(f, ::Type{T}, Pspace::S, Nspace::S,
-               [Espace::S], [Sspace::S], [Wspace::S]) where {T,S<:Union{Int,ElementarySpace}}
+               [Espace::S], [Sspace::S], [Wspace::S]) where {T, S<:Union{Int,ElementarySpace}}
                 
-Construct a PartitionFunctionTensor tensor based on the north, east, west and south spaces.
-The tensor elements are generated based on `f` and the element type is specified in `T`.
+Construct a `PartitionFunctionTensor` tensor based on the north, east, west and south spaces.
+The tensor elements are generated based on `f` and the element type and storage type is
+specified in `T`.
 """
 function PartitionFunctionTensor(
         f, ::Type{T},
@@ -80,13 +81,15 @@ const PEPSTensor{S <: ElementarySpace} = AbstractTensorMap{<:Any, S, 1, 4}
 
 """
     PEPSTensor(f, ::Type{T}, Pspace::S, Nspace::S,
-               [Espace::S], [Sspace::S], [Wspace::S]) where {T,S<:Union{Int,ElementarySpace}}
+               [Espace::S], [Sspace::S], [Wspace::S]) where {T, S<:Union{Int,ElementarySpace}}
                 
 Construct a PEPS tensor based on the physical, north, east, south and west spaces.
-The tensor elements are generated based on `f` and the element type is specified in `T`.
+The tensor elements are generated based on `f` and the element type and storage type
+is specified in `T`.
 """
 function PEPSTensor(
-        f, ::Type{T},
+        f,
+        ::Type{T},
         Pspace::S,
         Nspace::S, Espace::S = Nspace, Sspace::S = Nspace', Wspace::S = Espace',
     ) where {T, S <: ElementarySpace}
@@ -155,7 +158,7 @@ herm_depth(t::PEPOTensor) = permute(t', ((5, 6), (3, 2, 1, 4)))
 Fuse the physical indices of a PEPO tensor, obtaining a PEPS tensor.
 """
 function fuse_physicalspaces(O::PEPOTensor)
-    F = isomorphism(Int, fuse(codomain(O)), codomain(O))
+    F = isomorphism(TensorKit.similarstoragetype(O, Int), fuse(codomain(O)), codomain(O))
     return F * O, F
 end
 
diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl
index f902a3874..1f68a9161 100644
--- a/src/operators/infinitepepo.jl
+++ b/src/operators/infinitepepo.jl
@@ -121,7 +121,7 @@ function initializePEPS(
     end
     Nspaces = repeat([vspace], size(T, 1), size(T, 2))
     Espaces = repeat([vspace], size(T, 1), size(T, 2))
-    return InfinitePEPS(Pspaces, Nspaces, Espaces)
+    return InfinitePEPS(randn, storagetype(T), Pspaces, Nspaces, Espaces)
 end
 
 ## Unit cell interface
diff --git a/src/operators/transfermatrix.jl b/src/operators/transfermatrix.jl
index 85199a75f..fd94825e1 100644
--- a/src/operators/transfermatrix.jl
+++ b/src/operators/transfermatrix.jl
@@ -130,13 +130,13 @@ virtualspace(O::InfiniteTransferMatrix, i, dir) = virtualspace(O[i], dir)
 """
     initialize_mps(
         f=randn,
-        T=scalartype(O),
+        T=storagetype(O),
         O::Union{InfiniteTransferPEPS,InfiniteTransferPEPO},
         virtualspaces::AbstractArray{<:ElementarySpace,1}
     )
     initialize_mps(
         f=randn,
-        T=scalartype(O),
+        T=storagetype(O),
         O::Union{MultilineTransferPEPS,MultilineTransferPEPO},
         virtualspaces::AbstractArray{<:ElementarySpace,2}
     )
@@ -144,11 +144,11 @@ virtualspace(O::InfiniteTransferMatrix, i, dir) = virtualspace(O[i], dir)
 Inialize a boundary MPS for the transfer operator `O` by specifying an array of virtual
 spaces consistent with the unit cell.
 """
-function initialize_mps(O::Union{InfiniteTransferMatrix, MultilineTransferMatrix}, arg) # initialize(f=randn, T=scalartype(O), O, ...)
-    return initialize_mps(randn, scalartype(O), O, arg)
+function initialize_mps(O::Union{InfiniteTransferMatrix, MultilineTransferMatrix}, arg; kwargs...) # initialize(f=randn, T=storagetype(O), O, ...)
+    return initialize_mps(randn, storagetype(O), O, arg; kwargs...)
 end
 function initialize_mps(
-        f, T, O::InfiniteTransferMatrix, virtualspaces::AbstractArray{S, 1}
+        f, T, O::InfiniteTransferMatrix, virtualspaces::AbstractArray{S, 1}; kwargs...
     ) where {S}
     return InfiniteMPS(
         [
@@ -157,24 +157,24 @@ function initialize_mps(
                 virtualspaces[_prev(i, end)] * _elementwise_dual(north_virtualspace(O, i)),
                 virtualspaces[mod1(i, end)],
             ) for i in 1:length(O)
-        ]
+        ]; kwargs...
     )
 end
 function initialize_mps(
-        f, T, O::MultilineTransferMatrix, virtualspaces::AbstractArray{S, 2}
+        f, T, O::MultilineTransferMatrix, virtualspaces::AbstractArray{S, 2}; kwargs...
     ) where {S}
     mpss = map(1:size(O, 1)) do r
-        return initialize_mps(f, T, O[r], virtualspaces[r, :])
+        return initialize_mps(f, T, O[r], virtualspaces[r, :]; kwargs...)
     end
     return MPSKit.Multiline(mpss)
 end
 function initialize_mps(
-        f, T, O::MultilineTransferMatrix, virtualspaces::AbstractArray{S, 1}
+        f, T, O::MultilineTransferMatrix, virtualspaces::AbstractArray{S, 1}; kwargs...
     ) where {S}
-    return initialize_mps(f, T, O, repeat(virtualspaces, length(O), 1))
+    return initialize_mps(f, T, O, repeat(virtualspaces, length(O), 1); kwargs...)
 end
-function initialize_mps(f, T, O::MultilineTransferMatrix, V::ElementarySpace)
-    return initialize_mps(f, T, O, repeat([V], length(O), length(O[1])))
+function initialize_mps(f, T, O::MultilineTransferMatrix, V::ElementarySpace; kwargs...)
+    return initialize_mps(f, T, O, repeat([V], length(O), length(O[1])); kwargs...)
 end
 
 @doc """
diff --git a/src/states/infinitepeps.jl b/src/states/infinitepeps.jl
index 3d3a82f0f..d077f23e8 100644
--- a/src/states/infinitepeps.jl
+++ b/src/states/infinitepeps.jl
@@ -46,7 +46,7 @@ Create an `InfinitePEPS` by specifying the physical, north virtual and east virt
 of the PEPS tensor at each site in the unit cell as a matrix.
 """
 function InfinitePEPS(
-        f, T::Type{<:Number}, Pspaces::M, Nspaces::M, Espaces::M = Nspaces
+        f, T::Type, Pspaces::M, Nspaces::M, Espaces::M = Nspaces
     ) where {M <: AbstractMatrix{<:ElementarySpace}}
     size(Pspaces) == size(Nspaces) == size(Espaces) ||
         throw(ArgumentError("Input spaces should have equal sizes."))
@@ -108,7 +108,7 @@ end
 Create an InfinitePEPS by specifying its physical, north and east spaces and unit cell.
 """
 function InfinitePEPS(
-        f, T::Type{<:Number}, Pspace::S, vspaces...; unitcell::Tuple{Int, Int} = (1, 1)
+        f, T::Type, Pspace::S, vspaces...; unitcell::Tuple{Int, Int} = (1, 1)
     ) where {S <: ElementarySpace}
     return InfinitePEPS(
         f, T,
@@ -126,7 +126,7 @@ Base.eltype(::Type{InfinitePEPS{T}}) where {T} = T
 Base.eltype(A::InfinitePEPS) = eltype(typeof(A))
 
 Base.copy(A::InfinitePEPS) = InfinitePEPS(copy(unitcell(A)))
-function Base.similar(A::InfinitePEPS, T::Type{TorA} = scalartype(A)) where {TorA}
+function Base.similar(A::InfinitePEPS, T::Type = storagetype(A))
     return InfinitePEPS(map(t -> similar(t, T), unitcell(A)))
 end
 Base.repeat(A::InfinitePEPS, counts...) = InfinitePEPS(repeat(unitcell(A), counts...))