diff --git a/docs/src/cft.md b/docs/src/cft.md
index bde1504e..ac9c36a6 100644
--- a/docs/src/cft.md
+++ b/docs/src/cft.md
@@ -1,11 +1,12 @@
 # Conformal Field Theory Data
+
 TNRKit provides extensive tools for calculating conformal field theory data. Details about the implementation can be found in the TNRKit paper ([arxiv/2604.06922](https://arxiv.org/abs/2604.06922)).
 
-The core idea behind calculating the central charge, scaling dimensions, and conformal spins, is to calculate the spectrum of the fixed point tensor on a tube geometry. There are different ways to put fixed point tensors on a tube and the geometry of this tube is characterised by 3 parameters:
+The core idea behind calculating the central charge, scaling dimensions, and conformal spins is to calculate spectra of **transfer matrices** constructed from fixed point tensors on a tube geometry. There are different ways to put fixed point tensors on a tube, and the geometry of this tube is characterised by 3 parameters:
 $$[h, L, x]$$
-Where $h$ is the height of the tube, $L$ is the circumference, and $x$ is the horizontal shift. The higher the ratio $\frac{L}{h}$, the higher the resolution but also the more expensive the calculation.
+where $h$ is the height of the tube, $L$ is the circumference, and $x$ is the horizontal shift. The higher the ratio $\frac{L}{h}$, the higher the resolution but also the more expensive the calculation.
 
-To calculate cft data we provide the `CFTData` struct which can be used in the following ways:
+To calculate CFT data we provide the `CFTData` struct, which can be used in the following ways:
 
 ```julia
 CFTData(scheme; shape=[h, L, x])
@@ -15,18 +16,35 @@ CFTData(TA::TensorMap, TB::TensorMap; kwargs...) # 2x2 checkerboard unitcell
 
 The shapes we provide are: $[1, 1, 0]$, $[\sqrt{2}, 2\sqrt{2}, 0]$, $[1, 4, 1]$, $[1, 8, 1]$, $[\frac{4}{\sqrt{10}}, 2 \sqrt{10}, \frac{2}{\sqrt{10}}]$
 
-The last two of which require intermediate truncation steps, the parameters of which can be tuned by:
+The last two shapes require intermediate truncation steps, whose parameters can be tuned by:
+
 ```julia
 CFTData(scheme; shape=[1, 8, 1], trunc = trunc1, truncentanglement=trunc2)
 ```
 
-# CFTData struct
-The `CFTData` struct has two fields:
-- `central_charge`
-- `scaling_dimensions`
+## CFTData Struct
+
+The `CFTData` struct has three fields:
+
+- `central_charge`: a complex number which is the central charge of the CFT.
+- `modular_parameter`: a complex number which is the modular parameter of one tensor before building the tube transfer matrix.
+- `scaling_dimensions`: a `StructuredVector` storing the scaling dimensions and conformal spins of the primary fields and their descendants. It can be indexed like an `AbstractVector` (i.e. with scalars, slices, ...), or with sectors (e.g. `Z2Irrep(0)`), which will provide the scaling dimensions associated with that sector/charge.
+    To check which sectors you can index the `scaling_dimensions` with you can use `keys(scaling_dimensions)`.
+
+## How Fermionic CFT Data Is Extracted
 
-The `central_charge` can be either `missing` (when using the $[1, 1, 0]$ shape), or a number.
-The `scaling_dimensions` field is a `StructuredVector`.
+For a bosonic tensor network, `CFTData` builds a transfer matrix on the requested tube and diagonalizes it in each ordinary charge sector. For a fermionic tensor network, the same construction has one extra piece of global data: the **spin structure** around the tube. TNRKit extracts it by evaluating the same transfer matrix with two choices of boundary condition:
+
+- `:R`: periodic fermions around the tube, obtained with `pbc = true`. The macro `@tensor` automatically inserts a fermionic twist to take the *supertrace* across the tube.
+- `:NS`: antiperiodic fermions around the tube. Internally this is obtained by setting `pbc = false`, which explicitly inserts an additional fermionic twist to cancel the automatic supertrace twist, leaving the ordinary trace across the tube.
+
+The result is a `StructuredVector` whose keys are tuples `(spin_structure, charge)`. For fermionic systems without additional symmetries besides the fermion parity, there will be four sectors:
+
+```julia
+(:NS, FermionParity(0))
+(:NS, FermionParity(1))
+(:R,  FermionParity(0))
+(:R,  FermionParity(1))
+```
 
-The `scaling_dimensions` can be indexed like an `AbstractVector` (i.e. with scalars, slices, ...), or with sectors (e.g. `Z2Irrep(0)`), which will provide the scaling dimensions associated with that sector/charge.
-To check which sectors you can index the `scaling_dimensions` with you can use `keys(scaling_dimensions)`.
\ No newline at end of file
+The largest eigenvalue (corresponding to the identity field) should be in the NS even sector.
diff --git a/src/TNRKit.jl b/src/TNRKit.jl
index 7bc3a089..78f47824 100644
--- a/src/TNRKit.jl
+++ b/src/TNRKit.jl
@@ -15,6 +15,7 @@ using NonlinearSolve
 using Base.Threads
 using Combinatorics: permutations
 import TensorKitTensors.SpinOperators as SO
+import TensorKitTensors.FermionOperators as FO
 
 # stop criteria
 include("utility/stopping.jl")
@@ -125,7 +126,7 @@ export phi4_complex, phi4_complex_impϕ, phi4_complex_impϕdag, phi4_complex_imp
 
 include("models/quantum_1D.jl")
 export gate_to_tensor, vertical_stack_exp, vertical_stack_linear
-export quantum_ising_chain
+export quantum_ising_chain, kitaev_chain
 
 # utility functions
 include("utility/free_energy.jl")
diff --git a/src/models/quantum_1D.jl b/src/models/quantum_1D.jl
index fdce8ba1..97f69f4d 100644
--- a/src/models/quantum_1D.jl
+++ b/src/models/quantum_1D.jl
@@ -74,6 +74,11 @@ end
 # =================================================
 
 """
+    quantum_ising_chain(dt::Float64; kwargs...)
+    quantum_ising_chain(elt::Type{<:Number}, dt::Float64; kwargs...)
+    quantum_ising_chain(symm::Type{<:Sector}, dt::Float64; kwargs...)
+    quantum_ising_chain(elt::Type{<:Number}, symm::Type{<:Sector}, dt::Float64; J::Float64=1.0, g::Float64=0.0)
+
 Partition function tensor for 1D transverse field Ising chain
 ```
     H(PBC) = -J ∑_i (σz_i σz_{i+1} + g σx_i)
@@ -95,3 +100,59 @@ quantum_ising_chain(elt::Type{<:Number}, dt::Float64; kwargs...) =
     quantum_ising_chain(elt, Trivial, dt; kwargs...)
 quantum_ising_chain(symm::Type{<:Sector}, dt::Float64; kwargs...) =
     quantum_ising_chain(ComplexF64, symm, dt; kwargs...)
+quantum_ising_chain(dt::Float64; kwargs...) =
+    quantum_ising_chain(ComplexF64, Trivial, dt; kwargs...)
+
+"""
+    kitaev_chain(dt::Float64; kwargs...)
+    kitaev_chain(elt::Type{<:Number}, dt::Float64; kwargs...)
+    kitaev_chain(symm::Type{<:Sector}, dt::Float64; kwargs...)
+    kitaev_chain(elt::Type{<:Number}, symm::Type{<:Sector}, dt::Float64; t::Float64=1.0, Δ::Float64=1.0, V::Float64=0.0, µ::Float64=0.0)
+
+Partition function tensor for 1D Kitaev chain model
+```
+    H = ∑_i [
+        (-t) (c†_i c_{i+1} + h.c.) + Δ (c_i c_{i+1} + h.c.)
+        + V (n_i - 1/2) (n_{i+1} - 1/2) - μ(n_i - 1/2)
+    ]
+```
+It is related to the spin-1/2 Heisenberg XYZ model
+```
+    H = ∑_i (J_x Sx_i Sx_{i+1} + J_y Sy_i Sy_{i+1}
+            + J_z Sz_i Sz_{i+1} - h Sz_j)
+```
+by Jordan-Wigner transformation
+```
+    t = -(Jx + Jy) / 4,  V = Jz,
+    Δ = -(Jx - Jy) / 4,  μ = h.
+```
+Special Cases
+- t = 1, Δ = 1, V = 0, μ = 2:   transverse field Ising model
+- t = 1, Δ = 0, V ≠ 0, μ = 0:   Heisenberg XXZ model
+"""
+function kitaev_chain(
+        elt::Type{<:Number}, symm::Type{<:Sector}, dt::Float64;
+        t::Float64 = 1.0, Δ::Float64 = 1.0, V::Float64 = 0.0, µ::Float64 = 0.0
+    )
+    fpfm = FO.f_plus_f_min(elt, symm)
+    hopping = (-t) * (fpfm + fpfm')
+    num = FO.f_num(elt, symm)
+    unit = TensorKit.id(codomain(num, 1))
+    num = num - 0.5 * unit
+    interac = V * (num ⊗ num)
+    chempot = -µ * (num ⊗ unit + unit ⊗ num) / 2
+    gate = hopping + interac + chempot
+    if Δ != 0
+        fmfm = FO.f_min_f_min(elt, symm)
+        pairing = Δ * (fmfm + fmfm')
+        gate = gate + pairing
+    end
+    gate = exp(-dt * gate)
+    return gate_to_tensor(gate)
+end
+kitaev_chain(elt::Type{<:Number}, dt::Float64; kwargs...) =
+    kitaev_chain(elt, Trivial, dt; kwargs...)
+kitaev_chain(symm::Type{<:Sector}, dt::Float64; kwargs...) =
+    kitaev_chain(ComplexF64, symm, dt; kwargs...)
+kitaev_chain(dt::Float64; kwargs...) =
+    kitaev_chain(ComplexF64, Trivial, dt; kwargs...)
diff --git a/src/schemes/looptnr.jl b/src/schemes/looptnr.jl
index 8c5db9f1..8905cd94 100644
--- a/src/schemes/looptnr.jl
+++ b/src/schemes/looptnr.jl
@@ -9,12 +9,8 @@ Loop Optimization for Tensor Network Renormalization
     $(FUNCTIONNAME)(unitcell_2x2::Matrix{T})
 
 # Running the algorithm
-    run!(::LoopTNR, trunc::TruncationStrategy, criterion::stopcrit, parameters::LoopParameters, finalizer::Finalizer[,
-              entanglement_criterion::stopcrit, finalize_beginning=true, verbosity=1])
-    
-    run!(::LoopTNR, trscheme::TruncationStrategy, criterion::stopcrit, parameters::LoopParameters; kwargs...)
-
-    run!(::LoopTNR, trscheme::TruncationStrategy, criterion::stopcrit[finalize_beginning=true, verbosity=1])
+    run!(::LoopTNR, trunc::TruncationStrategy, criterion::stopcrit, [parameters::LoopParameters], [finalizer::Finalizer];
+        [entanglement_criterion::stopcrit, finalize_beginning=true, verbosity=1])
 
 # LoopParameters
 See also: [`LoopParameters`](@ref)
@@ -490,7 +486,7 @@ function step!(
 end
 
 function run!(
-        scheme::LoopTNR, trscheme::TruncationStrategy,
+        scheme::LoopTNR, trunc::TruncationStrategy,
         criterion::stopcrit, loop_condition::LoopParameters,
         finalizer::Finalizer{E};
         entanglement_criterion = default_entanglement_criterion,
@@ -510,7 +506,7 @@ function run!(
 
         t = @elapsed while crit
             @infov 2 "Step $(steps + 1), data[end]: $(!isempty(data) ? data[end] : "empty")"
-            step!(scheme, trscheme, entanglement_criterion, loop_condition, verbosity)
+            step!(scheme, trunc, entanglement_criterion, loop_condition, verbosity)
             push!(data, finalizer.f!(scheme))
 
             steps += 1
@@ -522,21 +518,23 @@ function run!(
     return data
 end
 
-function run!(scheme, trscheme, criterion, loop_condition; kwargs...)
-    return run!(scheme, trscheme, criterion, loop_condition, default_Finalizer; kwargs...)
-end
-
-function run!(
-        scheme::LoopTNR, trscheme::TruncationStrategy, criterion::stopcrit;
-        finalize_beginning = true, verbosity = 1
-    )
-    loop_condition = LoopParameters()
-    return run!(
-        scheme, trscheme, criterion, loop_condition;
-        finalize_beginning = finalize_beginning,
-        verbosity = verbosity
-    )
-end
+# use default finalizer
+run!(
+    scheme::LoopTNR, trunc::TruncationStrategy,
+    criterion::stopcrit, loop_condition::LoopParameters; kwargs...
+) = run!(scheme, trunc, criterion, loop_condition, default_Finalizer; kwargs...)
+
+# use default loop parameters
+run!(
+    scheme::LoopTNR, trunc::TruncationStrategy,
+    criterion::stopcrit, finalizer::Finalizer; kwargs...
+) = run!(scheme, trunc, criterion, LoopParameters(), finalizer; kwargs...)
+
+# use both default loop paramater and default finalizer
+run!(
+    scheme::LoopTNR, trunc::TruncationStrategy,
+    criterion::stopcrit; kwargs...
+) = run!(scheme, trunc, criterion, LoopParameters(), default_Finalizer; kwargs...)
 
 function Base.show(io::IO, scheme::LoopTNR)
     println(io, "LoopTNR - Loop Tensor Network Renormalization")
diff --git a/src/utility/cft.jl b/src/utility/cft.jl
index 2cc09fef..d603ced5 100644
--- a/src/utility/cft.jl
+++ b/src/utility/cft.jl
@@ -1,5 +1,5 @@
 """
-    struct CFTData{E, K, V, A <: AbstractVector{E}}
+    struct CFTData{E, K, A <: AbstractVector{E}}
 
 A struct to hold conformal data extracted from a TNR scheme.
 
@@ -11,16 +11,16 @@ A struct to hold conformal data extracted from a TNR scheme.
 # Fields
     - `central_charge::E`: The central charge of the CFT.
     - `modular_parameter::E`: The elementary modular parameter of a single tensor.
-    - `scaling_dimensions::StructuredVector{E, K, V, A}`: The scaling dimensions of the CFT, organized in a `StructuredVector` where the sectors correspond to different spin sectors (or other quantum numbers) and the data contains the scaling dimensions within those sectors
+    - `scaling_dimensions::StructuredVector{E, K, A}`: The scaling dimensions of the CFT, organized in a `StructuredVector` where the sectors correspond to different spin sectors (or other quantum numbers) and the data contains the scaling dimensions within those sectors
 
 """
-struct CFTData{E, K, V, A <: AbstractVector{E}}
+struct CFTData{E, K, A <: AbstractVector{E}}
     "Central charge of the CFT."
     central_charge::E
     "Elementary modular parameter for one tensor"
     modular_parameter::E
     "Scaling dimensions of the CFT."
-    scaling_dimensions::StructuredVector{E, K, V, A}
+    scaling_dimensions::StructuredVector{E, K, A}
 end
 
 function Base.show(io::IO, data::CFTData)
@@ -73,16 +73,7 @@ with `unitcell` copies of `T` concatenated horizontally.
 `τ0` is the modular parameter of a single `T`.
 """
 function _scaling_dimensions(T::TensorMap{E, S, 2, 2}, τ0::Number; unitcell = 1) where {E, S}
-    indices = [[i, -i, -(i + unitcell), i + 1] for i in 1:unitcell]
-    indices[end][4] = 1
-
-    T = ncon(fill(T, unitcell), indices)
-    # restore leg convention
-    outinds = Tuple(collect(1:unitcell))
-    ininds = Tuple(collect((unitcell + 1):(2unitcell)))
-    T = permute(T, (outinds, ininds))
-
-    sv = StructuredVector(eig_vals(T))
+    sv = _rowtm_eigvals(T, unitcell)
     sv = filter(x -> real(x) > 0 && abs(x) > 1.0e-12, sv)
     isempty(sv) && throw(ArgumentError("No valid eigenvalues found in transfer matrix spectrum."))
 
@@ -100,7 +91,8 @@ function area_term(
         TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}; is_real = true
     ) where {E, S}
     I = sectortype(TA)
-    λ = first(leading_eigenvalue(CFTTransferMatrix(TA, TB, [2, 2, 0]), one(I)))
+    pbc = (BraidingStyle(I) == Fermionic()) ? false : true
+    λ = first(leading_eigenvalue(CFTTransferMatrix(TA, TB, [2, 2, 0]), one(I); pbc, Nh = 1))
     return is_real ? real(λ) : λ
 end
 
@@ -120,10 +112,6 @@ function spec(
         τ0::Number; Nh = 25, trunc = notrunc(), truncentanglement = notrunc()
     ) where {E, S}
     I = sectortype(TA)
-    if BraidingStyle(I) != Bosonic()
-        throw(ArgumentError("Sectors with non-Bosonic charge $I has not been implemented"))
-    end
-
     tm = CFTTransferMatrix(TA, TB, shape; trunc, truncentanglement)
     τ = modular_parameter(tm, τ0)
 
@@ -131,7 +119,11 @@ function spec(
     eigs = leading_eigenvalue(tm; Nh)
 
     # central charge
-    λ0 = eigs[one(I)][1]
+    λ0 = if BraidingStyle(I) == Fermionic()
+        eigs[(:NS, one(I))][1]
+    else
+        eigs[one(I)][1]
+    end
     area = shape[1] * shape[2]
     central_charge = 6 / pi / (imag(τ) - imag(τ0) * area / 4) * log(λ0)
 
@@ -158,8 +150,10 @@ end
 sigmoid(x) = 1 / (1 + exp(-x))
 logit(p) = log(p / (1 - p))
 function _find_λ0(TA, TB, shape)
-    charge = one(sectortype(TA))
-    λs = leading_eigenvalue(CFTTransferMatrix(TA, TB, shape), charge; Nh = 1)
+    I = sectortype(TA)
+    charge = one(I)
+    pbc = (BraidingStyle(I) == Fermionic()) ? false : true
+    λs = leading_eigenvalue(CFTTransferMatrix(TA, TB, shape), charge; pbc, Nh = 1)
     return real(first(λs))
 end
 
diff --git a/src/utility/structuredvector.jl b/src/utility/structuredvector.jl
index ac9425c8..94047f85 100644
--- a/src/utility/structuredvector.jl
+++ b/src/utility/structuredvector.jl
@@ -1,9 +1,9 @@
 """
-    StructuredVector{E, K, V, A} <: AbstractVector{E}
+    StructuredVector{E, K, A} <: AbstractVector{E}
 
 A vector whose elements are partitioned into named sectors.  Internally, data
-is stored as a flat `AbstractVector{E}` and a `Dict{K, V}` maps each sector key
-to the indices that belong to it.
+is stored as a flat `AbstractVector{E}` and a `Dict{K, Vector{Int}}` maps each
+sector key to the indices that belong to it.
 
 Supports the `AbstractVector` interface (integer indexing, `length`, `eachindex`,
 …), sector-based access via `v[sector]`, and the full `Dict` key interface
@@ -14,15 +14,15 @@ broadcasting with scalars all preserve the sector structure.
 
     StructuredVector(sv::SectorVector)
     StructuredVector(dict::Dict{K, <:AbstractVector{E}}) where {K, E}
-    StructuredVector(data::AbstractVector{E}, structure::Dict{K, V})
+    StructuredVector(data::AbstractVector{E}, structure::Dict{K, Vector{Int}})
 
 - From a TensorKit `SectorVector`.
 - From a dictionary mapping sectors to their data vectors.
 - Directly from a flat data array and a sector‑index mapping.
 """
-struct StructuredVector{E, K, V, A <: AbstractVector{E}} <: AbstractVector{E}
+struct StructuredVector{E, K, A <: AbstractVector{E}} <: AbstractVector{E}
     data::A
-    structure::Dict{K, V}
+    structure::Dict{K, Vector{Int}}
 end
 
 function StructuredVector(sv::TensorKit.SectorVector)
@@ -84,6 +84,35 @@ Base.:/(x::Number, v::StructuredVector) = StructuredVector(x ./ v.data, v.struct
 
 Base.keys(v::StructuredVector) = keys(v.structure)
 
+function Base.vcat(v1::StructuredVector{<:Any, K1}, v2::StructuredVector{<:Any, K2}) where {K1, K2}
+    new_data = vcat(v1.data, v2.data)
+    n1 = length(v1)
+    K = promote_type(K1, K2)
+    new_structure = Dict{K, Vector{Int}}()
+    for (k, inds) in v1.structure
+        new_structure[k] = copy(inds)
+    end
+    for (k, inds) in v2.structure
+        adjusted = inds .+ n1
+        if haskey(new_structure, k)
+            append!(new_structure[k], adjusted)
+        else
+            new_structure[k] = adjusted
+        end
+    end
+    return StructuredVector(new_data, new_structure)
+end
+
+# Julia Base provides a specialized `reduce(::typeof(vcat), A)` that bypasses
+# pairwise `vcat` calls and uses `similar` + bulk copy internally.  That path
+# does not know about sector structure and would produce a plain `Vector`.
+# We override it to fall back to pair-wise reduction which preserves the
+# StructuredVector container.
+function Base.reduce(::typeof(vcat), A::AbstractVector{<:StructuredVector})
+    isempty(A) && throw(ArgumentError("reducing vcat over an empty collection is not supported"))
+    return foldl(vcat, A)
+end
+
 # -- Broadcasting support -----------------------------------------------------
 # Custom broadcast style so that element-wise operations preserve the
 # StructuredVector container (and therefore the sector structure).
@@ -111,3 +140,15 @@ function Base.similar(bc::Broadcast.Broadcasted{StructuredVectorStyle}, ::Type{E
     end
     return StructuredVector(similar(sv.data, ElType), copy(sv.structure))
 end
+
+"""
+    mapkeys(f, v::StructuredVector)
+
+Return a new `StructuredVector` with each `structure` key transformed
+by a function `f`. The underlying `data` is shared with `v`.
+```
+"""
+function mapkeys(f, v::StructuredVector)
+    new_structure = Dict(f(k) => copy(inds) for (k, inds) in v.structure)
+    return StructuredVector(v.data, new_structure)
+end
diff --git a/src/utility/transfer_matrix.jl b/src/utility/transfer_matrix.jl
index bf0fd758..a4244444 100644
--- a/src/utility/transfer_matrix.jl
+++ b/src/utility/transfer_matrix.jl
@@ -58,26 +58,26 @@ _has_shape(tm, shapes) = any(s -> tm.shape ≈ s, shapes)
 # ===========================================================================
 
 # Make the struct callable: (tm)(x) applies the transfer matrix
-function (tm::CFTTransferMatrix)(x)
-    return _TMaction(tm, x)
+function (tm::CFTTransferMatrix)(x; pbc::Bool = true)
+    return _TMaction(tm, x; pbc)
 end
 
 # Dispatch the action on the shape
-function _TMaction(tm::CFTTransferMatrix{E, S}, x) where {E, S}
+function _TMaction(tm::CFTTransferMatrix{E, S}, x; pbc::Bool = true) where {E, S}
     if tm.shape ≈ [1, 2, 1]
-        return _TMaction_1x2_NtoS(tm, x)
+        return _TMaction_1x2_NtoS(tm, x; pbc)
     elseif tm.shape ≈ [2, 2, 0]
-        return _TMaction_2x2_NtoS(tm, x)
+        return _TMaction_2x2_NtoS(tm, x; pbc)
     elseif tm.shape ≈ [sqrt(2) / 2, sqrt(2), sqrt(2) / 2]
-        return _TMaction_2x2_NEtoSW(tm, x)
+        return _TMaction_2x2_NEtoSW(tm, x; pbc)
     elseif tm.shape ≈ [sqrt(2), sqrt(2), 0]
-        return _TMaction_2x1_NEtoSW(tm, x)
+        return _TMaction_2x1_NEtoSW(tm, x; pbc)
     elseif tm.shape ≈ [1, 4, 1]
-        return _TMaction_1x4_twist(tm, x)
+        return _TMaction_1x4_twist(tm, x; pbc)
     elseif _has_shape(tm, _SHAPES_140)
-        return _TMaction_1x4(tm, x)
+        return _TMaction_1x4(tm, x; pbc)
     elseif _has_shape(tm, _SHAPES_2gates)
-        return _TMaction_2gates(tm, x)
+        return _TMaction_2gates(tm, x; pbc)
     else
         error("Unsupported transfer matrix shape: $(tm.shape).")
     end
@@ -93,9 +93,10 @@ Action of [1, 2, 1] transfer matrix.
     1'  2'
 ```
 """
-function _TMaction_1x2_NtoS(tm::CFTTransferMatrix{E, S}, x) where {E, S}
+function _TMaction_1x2_NtoS(tm::CFTTransferMatrix{E, S}, x; pbc::Bool = true) where {E, S}
+    TB′ = pbc ? tm.TB : twist(tm.TB, (2, 4))
     @tensor fx[-1 -2; -3] :=
-        tm.TA[a -2; 1 b] * twist(tm.TB, (2, 4))[b -1; 2 a] * x[1 2; -3]
+        tm.TA[a -2; 1 b] * TB′[b -1; 2 a] * x[1 2; -3]
     return fx
 end
 
@@ -111,10 +112,12 @@ Action of [2, 2, 0] transfer matrix.
         1'  2'
 ```
 """
-function _TMaction_2x2_NtoS(tm::CFTTransferMatrix{E, S}, x) where {E, S}
+function _TMaction_2x2_NtoS(tm::CFTTransferMatrix{E, S}, x; pbc::Bool = true) where {E, S}
+    TA′ = pbc ? tm.TA : twist(tm.TA, 1)
+    TB′ = pbc ? tm.TB : twist(tm.TB, 1)
     @tensor begin
-        fx[-1 -2; -3] := twist(tm.TA, 1)[c -1; 1 m] * x[1 2; -3] * tm.TB[m -2; 2 c]
-        fx[-1 -2; -3] := twist(tm.TB, 1)[c -1; 1 m] * fx[1 2; -3] * tm.TA[m -2; 2 c]
+        fx[-1 -2; -3] := TA′[c -1; 1 m] * x[1 2; -3] * tm.TB[m -2; 2 c]
+        fx[-1 -2; -3] := TB′[c -1; 1 m] * fx[1 2; -3] * tm.TA[m -2; 2 c]
     end
     return fx
 end
@@ -131,12 +134,14 @@ Action of [√2/2, √2, √2/2] transfer matrix.
         1'  2'
 ```
 """
-function _TMaction_2x2_NEtoSW(tm::CFTTransferMatrix{E, S}, x) where {E, S}
+function _TMaction_2x2_NEtoSW(tm::CFTTransferMatrix{E, S}, x; pbc::Bool = true) where {E, S}
+    TA′ = pbc ? tm.TA : twist(tm.TA, 2)
+    TB′ = pbc ? tm.TB : twist(tm.TB, 2)
     @tensor begin
         fx[-1 -2 -3 -4; -5] := tm.TB[-2 -3; a b] * x[-1 a b -4; -5]
         fx[-1 -2 -3 -4; -5] := tm.TA[-1 -2; a b] * fx[a b -3 -4; -5]
-        fx[-1 -2 -3 -4; -5] := twist(tm.TA, 2)[-3 -4; a b] * fx[-1 -2 a b; -5]
-        fx[-1 -2 -3 -4; -5] := twist(tm.TB, 2)[-4 -1; a b] * fx[-3 a b -2; -5]
+        fx[-1 -2 -3 -4; -5] := TA′[-3 -4; a b] * fx[-1 -2 a b; -5]
+        fx[-1 -2 -3 -4; -5] := TB′[-4 -1; a b] * fx[-3 a b -2; -5]
     end
     return fx
 end
@@ -153,10 +158,11 @@ Action of [√2, √2, 0] transfer matrix
         1'
 ```
 """
-function _TMaction_2x1_NEtoSW(tm::CFTTransferMatrix{E, S}, x) where {E, S}
+function _TMaction_2x1_NEtoSW(tm::CFTTransferMatrix{E, S}, x; pbc::Bool = true) where {E, S}
+    TB′ = pbc ? tm.TB : twist(tm.TB, (2, 4))
     @tensor begin
         fx[-1 -2; -3] := tm.TA[-1 -2; a b] * x[a b; -3]
-        fx[-1 -2; -3] := twist(tm.TB, (2, 4))[-2 -1; b a] * fx[a b; -3]
+        fx[-1 -2; -3] := TB′[-2 -1; b a] * fx[a b; -3]
     end
     return fx
 end
@@ -177,13 +183,14 @@ Action of [√2, 2√2, 0] transfer matrix
 First appeared in Chenfeng Bao's thesis: http://hdl.handle.net/10012/14674.
 """
 function _TMaction_2gates(
-        tm::CFTTransferMatrix{E, S}, x::TensorMap{E, S, 4, 1}
+        tm::CFTTransferMatrix{E, S}, x::TensorMap{E, S, 4, 1}; pbc::Bool = true
     ) where {E, S}
+    TA′ = pbc ? tm.TA : twist(tm.TA, (2, 4))
     @tensor begin
-        fx[-1 -2 -3 -4; -5] := tm.TB[-1 -2; 1 2] * x[1 2 3 4; -5] * tm.TB[-3 -4; 3 4]
-        fx[-1 -2 -3 -4; -5] := tm.TA[-3 -4; 2 3] * fx[1 2 3 4; -5] * tm.TA[-1 -2; 4 1]
+        fx[-1 -2 -3 -4; -5] := tm.TB[-1 -2; 1 2] * tm.TB[-3 -4; 3 4] * x[1 2 3 4; -5]
+        fx[-1 -2 -3 -4; -5] := tm.TA[-2 -3; 2 3] * TA′[-4 -1; 4 1] * fx[1 2 3 4; -5]
     end
-    return permute(fx, ((2, 3, 4, 1), (5,)))
+    return fx
 end
 
 """
@@ -197,11 +204,12 @@ Action of [1, 4, 0] transfer matrix. Only valid when TA = TB.
 ```
 """
 function _TMaction_1x4(
-        tm::CFTTransferMatrix{E, S}, x::TensorMap{E, S, 4, 1}
+        tm::CFTTransferMatrix{E, S}, x::TensorMap{E, S, 4, 1}; pbc::Bool = true
     ) where {E, S}
-    return @tensor TTTTx[-1 -2 -3 -4; -5] := x[1 2 3 4; -5] *
+    TB′ = pbc ? tm.TB : twist(tm.TB, 4)
+    return @tensor fx[-1 -2 -3 -4; -5] := x[1 2 3 4; -5] *
         tm.TA[41 -1; 1 12] * tm.TB[12 -2; 2 23] *
-        tm.TA[23 -3; 3 34] * tm.TB[34 -4; 4 41]
+        tm.TA[23 -3; 3 34] * TB′[34 -4; 4 41]
 end
 
 """
@@ -215,10 +223,13 @@ Action of [1, 4, 1] transfer matrix.
 ```
 """
 function _TMaction_1x4_twist(
-        tm::CFTTransferMatrix{E, S}, x::TensorMap{E, S, 4, 1}
+        tm::CFTTransferMatrix{E, S}, x::TensorMap{E, S, 4, 1}; pbc::Bool = true
     ) where {E, S}
-    TTTTx = _TMaction_1x4(tm, x)
-    return permute(TTTTx, ((2, 3, 4, 1), (5,)))
+    TB′ = pbc ? tm.TB : twist(tm.TB, (2, 4))
+    @tensor fx[-1 -2 -3 -4; -5] := x[1 2 3 4; -5] *
+        tm.TA[41 -2; 1 12] * tm.TB[12 -3; 2 23] *
+        tm.TA[23 -4; 3 34] * TB′[34 -1; 4 41]
+    return fx
 end
 
 # Utility to renormalize tensors for [1, 8, 1] and [4/√10, 2√10, 2/√10]
@@ -339,32 +350,81 @@ Returns a vector of eigenvalues.
 function leading_eigenvalue(tm::CFTTransferMatrix{E, S}; Nh::Int = 25) where {E, S}
     @assert Nh >= 1
     I = sectortype(tm.TA)
-    if BraidingStyle(I) != Bosonic()
-        throw(ArgumentError("Sectors with non-Bosonic charge $I has not been implemented"))
-    end
-
     xspace = domain(tm)
-    I = sectortype(tm.TA)
-    d = Dict{I, Vector{ComplexF64}}()
-    for charge in sectors(fuse(xspace))
-        vals = leading_eigenvalue(tm, charge; Nh = Nh)
-        isempty(vals) && continue
-        d[charge] = vals
+    if BraidingStyle(I) == Fermionic()
+        d = Dict{Tuple{Symbol, I}, Vector{ComplexF64}}()
+        for pbc in (false, true), charge in sectors(fuse(xspace))
+            vals = leading_eigenvalue(tm, charge; pbc, Nh)
+            isempty(vals) && continue
+            bc_label = pbc ? :R : :NS
+            d[(bc_label, charge)] = vals
+        end
+    else
+        d = Dict{I, Vector{ComplexF64}}()
+        for charge in sectors(fuse(xspace))
+            vals = leading_eigenvalue(tm, charge; pbc = true, Nh)
+            isempty(vals) && continue
+            d[charge] = vals
+        end
     end
     return StructuredVector(d)
 end
 
-function leading_eigenvalue(tm::CFTTransferMatrix{E, S}, charge; Nh::Int = 1) where {E, S}
+function leading_eigenvalue(
+        tm::CFTTransferMatrix{E, S}, charge; pbc::Bool = true, Nh::Int = 1
+    ) where {E, S}
     @assert Nh >= 1
     I = sectortype(tm.TA)
+    if !pbc
+        BraidingStyle(I) != Bosonic() || error("pbc = false only applies to fermionic tensors.")
+    end
     V = (I == Trivial) ? field(tm.TA)^1 : Vect[I](charge => 1)
     x = ones(domain(tm) ← V)
     dim(x) == 0 && error("$charge is not allowed by the transfer matrix.")
+    tm′(x) = tm(x; pbc)
     spec, _, info = eigsolve(
-        tm, x, Nh, :LM; krylovdim = 40, maxiter = 100, tol = 1.0e-12, verbosity = 0
+        tm′, x, Nh, :LM; krylovdim = 40, maxiter = 100, tol = 1.0e-12, verbosity = 0
     )
     if info.converged == 0
-        @warn "CFTTransferMatrix eigensolver did not converge in sector $charge."
+        @warn "CFTTransferMatrix eigensolver did not converge in sector $charge with pbc = $pbc."
     end
     return filter(x -> abs(real(x)) ≥ 1.0e-12, spec)
 end
+
+# ===========================================================================
+#  Special treatment: [1, L, 0] without `eigsolve`
+# ===========================================================================
+
+"""
+    _row_transfer_matrix(T::AbstractTensorMap, unitcell::Int; pbc::Bool = true)
+
+Build a row transfer matrix from `unitcell` copies of the tensor `T`.
+"""
+function _row_transfer_matrix(
+        T::AbstractTensorMap{E, S, 2, 2}, unitcell::Int; pbc::Bool = true
+    ) where {E, S}
+    indices = [[i, -i, -(i + unitcell), i + 1] for i in 1:unitcell]
+    indices[end][4] = 1
+    tensors = fill(T, unitcell)
+    if !pbc
+        tensors[end] = twist(T, 4)
+    end
+    Tcontracted = ncon(tensors, indices)
+    outinds = ntuple(i -> i, unitcell)
+    ininds = ntuple(i -> unitcell + i, unitcell)
+    return permute(Tcontracted, (outinds, ininds))
+end
+
+function _rowtm_eigvals(T::AbstractTensorMap{E, S, 2, 2}, unitcell::Int) where {E, S}
+    if BraidingStyle(sectortype(T)) == Fermionic()
+        # include contribution from both NS and R sectors
+        tm = twistdual(_row_transfer_matrix(T, unitcell; pbc = false), 1:unitcell)
+        ev_ns = mapkeys(k -> (:NS, k), StructuredVector(eig_vals(tm)))
+        tm = twistdual(_row_transfer_matrix(T, unitcell; pbc = true), 1:unitcell)
+        ev_rm = mapkeys(k -> (:R, k), StructuredVector(eig_vals(tm)))
+        return vcat(ev_ns, ev_rm)
+    else
+        tm = _row_transfer_matrix(T, unitcell; pbc = true)
+        return StructuredVector(eig_vals(tm))
+    end
+end
diff --git a/test/models/models.jl b/test/models/models.jl
index 07ef8878..d4ccbbd5 100644
--- a/test/models/models.jl
+++ b/test/models/models.jl
@@ -177,6 +177,14 @@ end
 end
 
 # (1 + 1)D quantum chains
+
+function cc_finalize!(scheme::LoopTNR)
+    n = finalize!(scheme)
+    τ, c = extract_tau_and_c(scheme.TA, scheme.TB)
+    @info "c = $c, τ = $τ"
+    return (n, τ, c)
+end
+
 @testset "Quantum Ising chain: $stack_alg stacking" for stack_alg in (:exponential, :linear)
     trunc_stack = truncerror(; rtol = 1.0e-9) & truncrank(16)
     if stack_alg == :exponential
@@ -191,9 +199,69 @@ end
         T = vertical_stack_linear(T, n, trunc_stack)
     end
     scheme = LoopTNR(T)
-    data = run!(scheme, truncrank(16), maxiter(16))
-    cft = CFTData(scheme)
-    central_charge = cft.central_charge
-    @test central_charge ≈ 0.5 atol = 1.0e-2
-    @info "Obtained central charge:\n$central_charge."
+    elt = scalartype(T)
+    finalizer = Finalizer(cc_finalize!, Tuple{elt, complex(elt), elt})
+    data = run!(scheme, truncrank(16), maxiter(16), finalizer; finalize_beginning = false)
+    @test last(data)[3] ≈ 0.5 atol = 1.0e-2
+end
+
+@testset "Quantum Kitaev chain: Ising case" begin
+    trunc_stack = truncerror(; rtol = 1.0e-9) & truncrank(16)
+    nfold = 7
+    dt = 1 / (2^nfold)
+    T = kitaev_chain(Float64, Trivial, dt; t = 1.0, Δ = 1.0, V = 0.0, μ = 2.0)
+    T = vertical_stack_exp(T, nfold, trunc_stack)
+
+    @info "CFT data from TRG"
+    scheme = TRG(T)
+    run!(scheme, truncrank(24), maxiter(8))
+    cft = CFTData(scheme; shape = [1, 1, 0])
+    c = cft.central_charge
+    sd = cft.scaling_dimensions
+    d_1 = real(sd[(:NS, FermionParity(0))][2])
+    d_f = real(sd[(:NS, FermionParity(1))][1])
+    d_e = real(sd[(:R, FermionParity(0))][1])
+    d_m = real(sd[(:R, FermionParity(1))][1])
+    @info "  shape [1, 1, 0]:\nΔ(1)=$d_1, Δ(f)=$d_f, Δ(e)=$d_e, Δ(m)=$d_m, c=$c"
+    @test d_1 ≈ 1 rtol = 2.0e-2
+    @test d_f ≈ 1 / 2 rtol = 2.0e-2
+    @test d_e ≈ 1 / 8 rtol = 2.0e-2
+    @test d_m ≈ 1 / 8 rtol = 2.0e-2
+    @test c ≈ 0.5 rtol = 2.0e-2
+
+    @info "CFT data from LoopTNR"
+    scheme = LoopTNR(T)
+    elt = scalartype(T)
+    finalizer = Finalizer(cc_finalize!, Tuple{elt, complex(elt), elt})
+    data = run!(scheme, truncrank(16), maxiter(8), finalizer; finalize_beginning = false)
+    @test last(data)[3] ≈ 0.5 atol = 1.0e-2
+    for shape in ([√2, 2√2, 0], [1, 4, 1], [1, 8, 1], [4 / √10, 2√10, 2 / √10])
+        cft = CFTData(scheme; shape = shape)
+        c = cft.central_charge
+        sd = cft.scaling_dimensions
+        d_1 = real(sd[(:NS, FermionParity(0))][2])
+        d_f = real(sd[(:NS, FermionParity(1))][1])
+        d_e = real(sd[(:R, FermionParity(0))][1])
+        d_m = real(sd[(:R, FermionParity(1))][1])
+        @info "  shape $shape:\nΔ(1)=$d_1, Δ(f)=$d_f, Δ(e)=$d_e, Δ(m)=$d_m, c=$c"
+        @test d_1 ≈ 1 rtol = 1.0e-2
+        @test d_f ≈ 1 / 2 rtol = 1.0e-2
+        @test d_e ≈ 1 / 8 rtol = 1.0e-2
+        @test d_m ≈ 1 / 8 rtol = 1.0e-2
+        @test c ≈ 0.5 rtol = 1.0e-2
+
+        shape[end] == 0 && continue
+        for sect in (
+                (:NS, FermionParity(0)), (:NS, FermionParity(1)),
+                (:R, FermionParity(0)), (:R, FermionParity(1)),
+            )
+            for v in sd[sect]
+                Δ, s = real(v), -imag(v)
+                Δ > 2.5 && break
+                s′ = (sect == (:NS, FermionParity(1))) ? s - 0.5 : s
+                @test isapprox(s′, round(s′); atol = 1.0e-2)
+            end
+            @info "  Conformal spins are integer (for 1, e, m) or half-integer (for f)."
+        end
+    end
 end