diff --git a/Project.toml b/Project.toml
index c721a55..20d53cf 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "DiagHamInterface"
 uuid = "2f27495d-c82b-46d5-b81b-1bb9aa58c4d5"
-version = "1.1.1"
+version = "1.2.0"
 authors = ["Andreas Feuerpfeil <development@manybodylab.com>"]
 
 [deps]
diff --git a/docs/Project.toml b/docs/Project.toml
index a78e2d9..e3ac7a3 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -6,4 +6,4 @@ DiagHamInterface = "2f27495d-c82b-46d5-b81b-1bb9aa58c4d5"
 [compat]
 Documenter = "1"
 Literate = "2"
-DiagHamInterface = "1.0,1.1"
+DiagHamInterface = "1.0,1.1,1.2"
diff --git a/src/DiagHamInterface.jl b/src/DiagHamInterface.jl
index 33ab05f..d03c099 100644
--- a/src/DiagHamInterface.jl
+++ b/src/DiagHamInterface.jl
@@ -11,6 +11,11 @@ export read_matrix_from_txt
 export write_to_txt
 export write_matrix_elements
 
+export read_one_body_density
+export read_two_body_density
+export two_body_normal_to_density_density
+export density_density_to_two_body_normal
+
 using DelimitedFiles
 using SparseArrays
 using Preferences
@@ -26,5 +31,6 @@ include("utility/standards.jl")
 include("MatrixElements/read.jl")
 include("MatrixElements/write.jl")
 include("Hamiltonian/hamiltonian_txt.jl")
+include("FTIDensity/read_density.jl")
 
 end
diff --git a/src/FTIDensity/read_density.jl b/src/FTIDensity/read_density.jl
new file mode 100644
index 0000000..f91956f
--- /dev/null
+++ b/src/FTIDensity/read_density.jl
@@ -0,0 +1,265 @@
+function _parse_density_header(filename::AbstractString)
+    header_line = open(filename) do io
+        for line in eachline(io)
+            startswith(strip(line), '#') && return line
+        end
+        error("No header line found in $filename")
+    end
+
+    is_two_body = occursin("c^+ c^+ c c", header_line)
+
+    tokens = split(header_line)
+    # Remove leading '#' and drop everything from the first '<…>' group onward
+    filter!(t -> t != "#", tokens)
+    bracket_start = findfirst(t -> startswith(t, "<"), tokens)
+    if !isnothing(bracket_start)
+        resize!(tokens, bracket_start - 1)
+    end
+    Ncols_index = length(tokens)
+
+    has_spin = any(t -> t == "spin" || startswith(t, "spin"), tokens)
+
+    # Catch unsupported modes
+    if any(occursin(t, "psi_i") || occursin(t, "phi_j") for t in tokens)
+        error("Off-diagonal mode (--off-diagonal) output is not supported")
+    end
+    if any(t -> t == "kz" || startswith(t, "kz"), tokens)
+        error("3D lattice output is not supported")
+    end
+
+    return is_two_body, has_spin, Ncols_index
+end
+
+function _read_density_rows(filename::AbstractString, Ncols_index::Int)
+    index_rows = Vector{Int}[]
+    value_strs = String[]
+
+    open(filename) do io
+        first_header_seen = false
+        for line in eachline(io)
+            stripped = strip(line)
+            isempty(stripped) && continue
+            if startswith(stripped, '#')
+                first_header_seen = true
+                continue
+            end
+            first_header_seen || continue  # skip anything before first header
+            tokens = split(stripped)
+            length(tokens) < Ncols_index + 1 && continue
+            push!(index_rows, parse.(Int, tokens[1:Ncols_index]))
+            push!(value_strs, join(tokens[(Ncols_index + 1):end], " "))
+        end
+    end
+
+    isempty(index_rows) && return Matrix{Int}(undef, 0, Ncols_index), Number[]
+    index_matrix = reduce(hcat, index_rows)'  # (nrows, Ncols_index)
+    values = read_number.(value_strs)
+    return index_matrix, values
+end
+
+function _determine_type(values::AbstractVector, atol::Real)
+    isempty(values) && return Float64
+    return any(v -> abs(imag(v)) > atol, values) ? ComplexF64 : Float64
+end
+
+"""
+    read_one_body_density(filename; atol=1e-14) -> Array{T,6}
+
+Read a one-body density matrix file produced by DiagHam's FTIDensity program.
+
+Returns an array of shape `(Nk, Nband, Nspin, Nk, Nband, Nspin)` where
+`result[k1, b1, s1, k2, b2, s2] = ⟨c†_{k1-1, b1-1, s1-1} c_{k2-1, b2-1, s2-1}⟩`.
+
+The momentum index `k` is linearized as `k = kx*Nky + ky` (0-based in file,
+1-based in the returned array). Nkx, Nky, Nband, Nspin are inferred from the
+data. For single-band files Nband=1 and Nspin=1. For files with a `spin` column
+(NbrBands=4 or 6) Nspin=2.
+
+Since FTIDensity only writes diagonal-in-k entries, all off-diagonal-k entries
+are zero.
+"""
+function read_one_body_density(filename::AbstractString; atol::Real = 1.0e-14)
+    is_two_body, has_spin, Ncols_index = _parse_density_header(filename)
+    is_two_body && error("File contains two-body density; use read_two_body_density instead")
+
+    index_matrix, values = _read_density_rows(filename, Ncols_index)
+
+    if isempty(values)
+        T = Float64
+        return zeros(T, 0, 1, 1, 0, 1, 1)
+    end
+
+    T = _determine_type(values, atol)
+
+    kx_col = index_matrix[:, 1]
+    ky_col = index_matrix[:, 2]
+    Nkx = maximum(kx_col) + 1
+    Nky = maximum(ky_col) + 1
+    Nk = Nkx * Nky
+
+    if Ncols_index == 2
+        # 1-band: kx ky val
+        Nband = 1
+        Nspin = 1
+        result = zeros(T, Nk, Nband, Nspin, Nk, Nband, Nspin)
+        for (i, val) in enumerate(values)
+            k = kx_col[i] * Nky + ky_col[i] + 1
+            v = T == Float64 ? real(val) : convert(T, val)
+            result[k, 1, 1, k, 1, 1] = v
+        end
+    elseif !has_spin
+        # 2/3-band: kx ky sigma_cr sigma_an val
+        sigcr_col = index_matrix[:, 3]
+        sigan_col = index_matrix[:, 4]
+        Nband = maximum(max.(sigcr_col, sigan_col)) + 1
+        Nspin = 1
+        result = zeros(T, Nk, Nband, Nspin, Nk, Nband, Nspin)
+        for (i, val) in enumerate(values)
+            k = kx_col[i] * Nky + ky_col[i] + 1
+            bcr = sigcr_col[i] + 1
+            ban = sigan_col[i] + 1
+            v = T == Float64 ? real(val) : convert(T, val)
+            result[k, bcr, 1, k, ban, 1] = v
+        end
+    else
+        # 4/6-band: kx ky spin sigma_cr sigma_an val
+        spin_col = index_matrix[:, 3]
+        sigcr_col = index_matrix[:, 4]
+        sigan_col = index_matrix[:, 5]
+        Nband = maximum(max.(sigcr_col, sigan_col)) + 1
+        Nspin = maximum(spin_col) + 1
+        result = zeros(T, Nk, Nband, Nspin, Nk, Nband, Nspin)
+        for (i, val) in enumerate(values)
+            k = kx_col[i] * Nky + ky_col[i] + 1
+            s = spin_col[i] + 1
+            bcr = sigcr_col[i] + 1
+            ban = sigan_col[i] + 1
+            v = T == Float64 ? real(val) : convert(T, val)
+            result[k, bcr, s, k, ban, s] = v
+        end
+    end
+
+    return result
+end
+
+"""
+    read_two_body_density(filename; atol=1e-14) -> Array{T,12}
+
+Read a two-body density matrix file produced by DiagHam's FTIDensity program
+(output of the `--rhorho` flag).
+
+Returns an array of shape `(Nk,Nb,Ns, Nk,Nb,Ns, Nk,Nb,Ns, Nk,Nb,Ns)` where
+`result[k1,b1,s1, k2,b2,s2, k3,b3,s3, k4,b4,s4] = ⟨c†_1 c†_2 c_3 c_4⟩`.
+All indices are 1-based; k is linearized as `k = kx*Nky + ky + 1`.
+"""
+function read_two_body_density(filename::AbstractString; atol::Real = 1.0e-14)
+    is_two_body, has_spin, Ncols_index = _parse_density_header(filename)
+    !is_two_body && error("File contains one-body density; use read_one_body_density instead")
+
+    index_matrix, values = _read_density_rows(filename, Ncols_index)
+
+    if isempty(values)
+        T = Float64
+        return zeros(T, ntuple(_ -> 1, 12)...)
+    end
+
+    T = _determine_type(values, atol)
+
+    stride = has_spin ? 4 : 3  # columns per particle: kx, ky, [spin,] sigma
+
+    kx_cols = [index_matrix[:, 1 + (p - 1) * stride] for p in 1:4]
+    ky_cols = [index_matrix[:, 2 + (p - 1) * stride] for p in 1:4]
+    if has_spin
+        spin_cols = [index_matrix[:, 3 + (p - 1) * stride] for p in 1:4]
+        sigma_cols = [index_matrix[:, 4 + (p - 1) * stride] for p in 1:4]
+    else
+        sigma_cols = [index_matrix[:, 3 + (p - 1) * stride] for p in 1:4]
+    end
+
+    Nkx = maximum(maximum.(kx_cols)) + 1
+    Nky = maximum(maximum.(ky_cols)) + 1
+    Nk = Nkx * Nky
+    Nband = maximum(maximum.(sigma_cols)) + 1
+    Nspin = has_spin ? maximum(maximum.(spin_cols)) + 1 : 1
+
+    result = zeros(T, Nk, Nband, Nspin, Nk, Nband, Nspin, Nk, Nband, Nspin, Nk, Nband, Nspin)
+
+    for (i, val) in enumerate(values)
+        ks = [kx_cols[p][i] * Nky + ky_cols[p][i] + 1 for p in 1:4]
+        bs = [sigma_cols[p][i] + 1 for p in 1:4]
+        ss = has_spin ? [spin_cols[p][i] + 1 for p in 1:4] : [1, 1, 1, 1]
+        v = T == Float64 ? real(val) : convert(T, val)
+        result[ks[1], bs[1], ss[1], ks[2], bs[2], ss[2], ks[3], bs[3], ss[3], ks[4], bs[4], ss[4]] = v
+    end
+
+    return result
+end
+
+"""
+    two_body_normal_to_density_density(rho2, rho1; statistics=:fermi) -> Array{T,12}
+
+Convert the normal-ordered two-body density matrix `⟨c†_1 c†_2 c_3 c_4⟩` to
+the density-density ordering `⟨c†_1 c_3 c†_2 c_4⟩` using commutation relations.
+
+For fermions (`statistics=:fermi`, anticommutator `{c,c†}=δ`):
+    ⟨c†_1 c_3 c†_2 c_4⟩ = δ_{23} ⟨c†_1 c_4⟩ - ⟨c†_1 c†_2 c_3 c_4⟩
+
+For bosons (`statistics=:bose`, commutator `[b,b†]=δ`):
+    ⟨b†_1 b_3 b†_2 b_4⟩ = δ_{23} ⟨b†_1 b_4⟩ + ⟨b†_1 b†_2 b_3 b_4⟩
+
+`rho2` has shape `(Nk,Nb,Ns, Nk,Nb,Ns, Nk,Nb,Ns, Nk,Nb,Ns)` (from
+`read_two_body_density`). `rho1` has shape `(Nk,Nb,Ns, Nk,Nb,Ns)` (from
+`read_one_body_density`). The returned array has the same shape as `rho2`, with
+index semantics `[k1,b1,s1, k3,b3,s3, k2,b2,s2, k4,b4,s4]`.
+"""
+function two_body_normal_to_density_density(
+        rho2::Array{T, 12}, rho1::Array{T, 6}; statistics::Symbol = :fermi
+    ) where {T}
+    coeff = _rho2_permuted_coeff(statistics)
+    # Permute rho2: swap particle-2 block (dims 4-6) and particle-3 block (dims 7-9)
+    rho_dd = coeff .* permutedims(rho2, [1, 2, 3, 7, 8, 9, 4, 5, 6, 10, 11, 12])
+    _add_delta_rho1!(rho_dd, rho1)
+    return rho_dd
+end
+
+"""
+    density_density_to_two_body_normal(rho_dd, rho1; statistics=:fermi) -> Array{T,12}
+
+Convert the density-density ordering `⟨c†_1 c_3 c†_2 c_4⟩` back to the
+normal-ordered two-body density matrix `⟨c†_1 c†_2 c_3 c_4⟩`.
+
+For fermions: `⟨c†_1 c†_2 c_3 c_4⟩ = δ_{23} ⟨c†_1 c_4⟩ - ⟨c†_1 c_3 c†_2 c_4⟩`
+For bosons:   `⟨b†_1 b†_2 b_3 b_4⟩ = δ_{23} ⟨b†_1 b_4⟩ + ⟨b†_1 b_3 b†_2 b_4⟩`
+
+`rho_dd` has index semantics `[k1,b1,s1, k3,b3,s3, k2,b2,s2, k4,b4,s4]` as
+returned by `two_body_normal_to_density_density`. The returned array has the
+same shape with index semantics `[k1,b1,s1, k2,b2,s2, k3,b3,s3, k4,b4,s4]`.
+"""
+function density_density_to_two_body_normal(
+        rho_dd::Array{T, 12}, rho1::Array{T, 6}; statistics::Symbol = :fermi
+    ) where {T}
+    coeff = _rho2_permuted_coeff(statistics)
+    # Same permutation: swap particle-2 block (dims 7-9 in rho_dd) and particle-3 block (dims 4-6)
+    rho2 = coeff .* permutedims(rho_dd, [1, 2, 3, 7, 8, 9, 4, 5, 6, 10, 11, 12])
+    _add_delta_rho1!(rho2, rho1)
+    return rho2
+end
+
+# Coefficient applied to the permuted rho2 (or rho_dd) in the conversion formula:
+#   fermi: rho_dd = -rho2_perm + δ·rho1  → coeff = -1
+#   bose:  rho_dd = +rho2_perm + δ·rho1  → coeff = +1
+function _rho2_permuted_coeff(statistics::Symbol)
+    statistics == :fermi && return -1
+    statistics == :bose  && return +1
+    error("Unknown statistics $statistics; expected :fermi or :bose")
+end
+
+# Add δ_{23} ⟨c†_1 c_4⟩ to result[k1,b1,s1, k2,b2,s2, k2,b2,s2, k4,b4,s4]
+# (dims 4-6 and 7-9 are the two α2=α3 blocks, dims 1-3=α1, 10-12=α4)
+function _add_delta_rho1!(result::Array{T, 12}, rho1::Array{T, 6}) where {T}
+    Nk, Nb, Ns = size(rho1, 1), size(rho1, 2), size(rho1, 3)
+    for s2 in 1:Ns, b2 in 1:Nb, k2 in 1:Nk
+        @views result[:, :, :, k2, b2, s2, k2, b2, s2, :, :, :] .+= rho1[:, :, :, :, :, :]
+    end
+    return result
+end
diff --git a/src/MatrixElements/read.jl b/src/MatrixElements/read.jl
index c0ba567..14252a2 100644
--- a/src/MatrixElements/read.jl
+++ b/src/MatrixElements/read.jl
@@ -57,8 +57,12 @@ function reinstate_indices(header, indices, coeffs) ## This function takes care
         # positions of headers with suffix s
         pos = findall(x -> suffixes[x] == s, 1:length(header))
         if !isempty(pos)
-            # local order within the group: sort by label string (stable)
-            local_order = sortperm(header[pos])
+            label_priority = h -> begin
+                base = split(h, "_")[1]
+                startswith(base, standard_band_label()) ? (0, base) :
+                    startswith(base, standard_valley_label()) ? (1, base) : (2, base)
+            end
+            local_order = sortperm(header[pos], by = label_priority)
             append!(perm, pos[local_order])
         end
     end
diff --git a/src/MatrixElements/write.jl b/src/MatrixElements/write.jl
index 91c5ac8..2eab881 100644
--- a/src/MatrixElements/write.jl
+++ b/src/MatrixElements/write.jl
@@ -74,6 +74,22 @@ function format_kinetic(label, indices, coeffs; full_single_particle::Bool = fal
     return label, indices, coeffs
 end
 
+function reorder_label_columns(label, indices)
+    label_priority = h -> begin
+        base = split(h, "_")[1]
+        startswith(base, standard_band_label()) ? (0, base) :
+            startswith(base, standard_valley_label()) ? (1, base) : (2, base)
+    end
+    suffixes = [parse(Int, split(h, "_")[2]) for h in label]
+    perm = Int[]
+    for s in unique(sort(suffixes))
+        pos = findall(x -> suffixes[x] == s, 1:length(label))
+        local_order = sortperm(label[pos], by = label_priority)
+        append!(perm, pos[local_order])
+    end
+    return label[perm], indices[:, perm]
+end
+
 """
     write_matrix_elements(label, indices, coeffs, file_name; dropband=false, full_single_particle=false)
 
@@ -88,6 +104,9 @@ function write_matrix_elements(label, indices, coeffs, file_name::String; dropba
     label, indices, coeffs = drop_bands(label, indices, coeffs; dropband = dropband)
     label, indices, coeffs = relabel_bands(label, indices, coeffs; full_single_particle = full_single_particle)
 
+    if N_body_interaction > 1
+        label, indices = reorder_label_columns(label, indices)
+    end
 
     # Finally store
     fid = open(file_name, "w")
diff --git a/src/utility/numbers.jl b/src/utility/numbers.jl
index 4907cff..5a10335 100644
--- a/src/utility/numbers.jl
+++ b/src/utility/numbers.jl
@@ -1,11 +1,11 @@
 """
-    format_with_precision(x; atol=0.0, maxdigits=typemax(Int))
+    format_with_precision(x; atol=eps(float(T)), maxdigits=typemax(Int))
 
 Format number `x` as a string with absolute precision `atol` and a maximum of 
 `maxdigits` significant digits. 
 Uses scientific notation for very small or large numbers.
 """
-function format_with_precision(x::T; atol = 0.0, maxdigits::Int = typemax(Int)) where {T <: Real}
+function format_with_precision(x::T; atol = eps(float(T)), maxdigits::Int = typemax(Int)) where {T <: Real}
     iszero(x) && return "0.0"
 
     if iszero(atol)
diff --git a/test/FTIDensity/test_read_density.jl b/test/FTIDensity/test_read_density.jl
new file mode 100644
index 0000000..8b14f9a
--- /dev/null
+++ b/test/FTIDensity/test_read_density.jl
@@ -0,0 +1,249 @@
+using DiagHamInterface
+using DiagHamInterface: _parse_density_header
+using Test
+
+# ── helpers ──────────────────────────────────────────────────────────────────
+
+function write_tmp(content::String)
+    path, io = mktemp()
+    write(io, content)
+    close(io)
+    return path
+end
+
+# ── header parsing ────────────────────────────────────────────────────────────
+
+@testset "Header parsing" begin
+    @testset "1-band" begin
+        f = write_tmp("# kx ky <c^+ c>\n0 0 1.0\n")
+        is2, hspin, ncols = _parse_density_header(f)
+        @test !is2
+        @test !hspin
+        @test ncols == 2
+    end
+
+    @testset "2-band" begin
+        f = write_tmp("# kx ky sigma sigma' <c^+ c>\n0 0 0 0 1.0\n")
+        is2, hspin, ncols = _parse_density_header(f)
+        @test !is2
+        @test !hspin
+        @test ncols == 4
+    end
+
+    @testset "4-band (spin)" begin
+        f = write_tmp("# kx ky spin sigma sigma' <c^+ c>\n0 0 0 0 0 1.0\n")
+        is2, hspin, ncols = _parse_density_header(f)
+        @test !is2
+        @test hspin
+        @test ncols == 5
+    end
+
+    @testset "two-body 2-band" begin
+        f = write_tmp("# kx1 ky1 sigma1 kx2 ky2 sigma2 kx3 ky3 sigma3 kx4 ky4 sigma4 <c^+ c^+ c c>\n0 0 0 0 0 0 0 0 0 0 0 0 1.0\n")
+        is2, hspin, ncols = _parse_density_header(f)
+        @test is2
+        @test !hspin
+        @test ncols == 12
+    end
+end
+
+# ── one-body reader ───────────────────────────────────────────────────────────
+
+@testset "read_one_body_density" begin
+    @testset "1-band, 2x2 lattice, real" begin
+        content = """# kx ky <c^+ c>
+        0 0 0.5
+        0 1 0.25
+        1 0 0.125
+        1 1 0.0625
+        # partial density sigma=0 = 0.9375
+        # total density = 0.9375
+        """
+        f = write_tmp(content)
+        rho = read_one_body_density(f)
+
+        @test eltype(rho) == Float64
+        @test size(rho) == (4, 1, 1, 4, 1, 1)
+
+        # k = kx*Nky + ky + 1  (Nky=2)
+        @test rho[1, 1, 1, 1, 1, 1] ≈ 0.5       # k=0*2+0+1=1
+        @test rho[2, 1, 1, 2, 1, 1] ≈ 0.25      # k=0*2+1+1=2
+        @test rho[3, 1, 1, 3, 1, 1] ≈ 0.125     # k=1*2+0+1=3
+        @test rho[4, 1, 1, 4, 1, 1] ≈ 0.0625    # k=1*2+1+1=4
+
+        # Off-diagonal k entries must be zero
+        @test rho[1, 1, 1, 2, 1, 1] == 0
+        @test rho[2, 1, 1, 3, 1, 1] == 0
+    end
+
+    @testset "2-band, 1x1 lattice, complex" begin
+        content = """# kx ky sigma sigma' <c^+ c>
+        0 0 0 0 1.0
+        0 0 0 1 (0.1,0.2)
+        0 0 1 0 (0.1,-0.2)
+        0 0 1 1 0.3
+        """
+        f = write_tmp(content)
+        rho = read_one_body_density(f)
+
+        @test eltype(rho) == ComplexF64
+        @test size(rho) == (1, 2, 1, 1, 2, 1)
+
+        @test rho[1, 1, 1, 1, 1, 1] ≈ 1.0
+        @test rho[1, 1, 1, 1, 2, 1] ≈ 0.1 + 0.2im
+        @test rho[1, 2, 1, 1, 1, 1] ≈ 0.1 - 0.2im
+        @test rho[1, 2, 1, 1, 2, 1] ≈ 0.3
+
+        # Hermiticity
+        @test rho[1, 1, 1, 1, 2, 1] ≈ conj(rho[1, 2, 1, 1, 1, 1])
+    end
+
+    @testset "4-band (spin), 1x1 lattice" begin
+        content = """# kx ky spin sigma sigma' <c^+ c>
+        0 0 0 0 0 0.8
+        0 0 0 0 1 0.0
+        0 0 0 1 0 0.0
+        0 0 0 1 1 0.2
+        0 0 1 0 0 0.6
+        0 0 1 0 1 0.0
+        0 0 1 1 0 0.0
+        0 0 1 1 1 0.4
+        """
+        f = write_tmp(content)
+        rho = read_one_body_density(f)
+
+        @test eltype(rho) == Float64
+        @test size(rho) == (1, 2, 2, 1, 2, 2)
+
+        @test rho[1, 1, 1, 1, 1, 1] ≈ 0.8   # spin=0, sigma=0→0
+        @test rho[1, 2, 1, 1, 2, 1] ≈ 0.2   # spin=0, sigma=1→1
+        @test rho[1, 1, 2, 1, 1, 2] ≈ 0.6   # spin=1, sigma=0→0
+        @test rho[1, 2, 2, 1, 2, 2] ≈ 0.4   # spin=1, sigma=1→1
+
+        # Cross-spin entries must be zero
+        @test rho[1, 1, 1, 1, 1, 2] == 0
+        @test rho[1, 1, 2, 1, 1, 1] == 0
+    end
+
+    @testset "Comment lines interspersed with data" begin
+        # DiagHam inserts partial-density comments between data blocks
+        content = """# kx ky sigma sigma' <c^+ c>
+        0 0 0 0 1.0
+        # partial density sigma=0 = (1.0,0.0)
+        0 0 1 1 0.5
+        # partial density sigma=1 = (0.5,0.0)
+        # total density = (1.5,0.0)
+        """
+        f = write_tmp(content)
+        rho = read_one_body_density(f)
+        @test size(rho) == (1, 2, 1, 1, 2, 1)
+        @test rho[1, 1, 1, 1, 1, 1] ≈ 1.0
+        @test rho[1, 2, 1, 1, 2, 1] ≈ 0.5
+    end
+
+    @testset "Real detection: purely real values → Float64" begin
+        content = "# kx ky <c^+ c>\n0 0 0.7\n"
+        f = write_tmp(content)
+        rho = read_one_body_density(f)
+        @test eltype(rho) == Float64
+    end
+
+    @testset "Complex detection: nonzero imaginary → ComplexF64" begin
+        content = "# kx ky <c^+ c>\n0 0 (0.7,0.1)\n"
+        f = write_tmp(content)
+        rho = read_one_body_density(f)
+        @test eltype(rho) == ComplexF64
+    end
+end
+
+# ── two-body reader ───────────────────────────────────────────────────────────
+
+@testset "read_two_body_density" begin
+    @testset "2-band, 2x1 lattice" begin
+        # kx1 ky1 sig1 kx2 ky2 sig2 kx3 ky3 sig3 kx4 ky4 sig4 val
+        # Momentum conservation: kx4 = (kx1+kx2-kx3) mod Nkx
+        content = """# kx1 ky1 sigma1 kx2 ky2 sigma2 kx3 ky3 sigma3 kx4 ky4 sigma4 <c^+ c^+ c c>
+        0 0 0 1 0 1 0 0 0 1 0 1 0.5
+        1 0 0 0 0 1 1 0 0 0 0 1 0.25
+        """
+        f = write_tmp(content)
+        rho2 = read_two_body_density(f)
+
+        # Nkx=2, Nky=1 → Nk=2; Nband=2; Nspin=1
+        @test size(rho2) == (2, 2, 1, 2, 2, 1, 2, 2, 1, 2, 2, 1)
+        @test eltype(rho2) == Float64
+
+        # k=kx*Nky+ky+1; Nky=1
+        # row 1: k1=1,b1=1,s1=1; k2=2,b2=2,s2=1; k3=1,b3=1,s3=1; k4=2,b4=2,s4=1
+        @test rho2[1, 1, 1, 2, 2, 1, 1, 1, 1, 2, 2, 1] ≈ 0.5
+        # row 2: k1=2,b1=1,s1=1; k2=1,b2=2,s2=1; k3=2,b3=1,s3=1; k4=1,b4=2,s4=1
+        @test rho2[2, 1, 1, 1, 2, 1, 2, 1, 1, 1, 2, 1] ≈ 0.25
+    end
+end
+
+# ── commutation converters ────────────────────────────────────────────────────
+
+@testset "Commutation converters" begin
+    # Product state |ψ⟩ = c†_a c†_b |0⟩, a=(k=1,b=1,s=1), b=(k=2,b=1,s=1)
+    # One-body: ⟨c†_a c_a⟩=1, ⟨c†_b c_b⟩=1, others=0
+    # Two-body: ⟨c†_a c†_b c_b c_a⟩=1, ⟨c†_a c†_b c_a c_b⟩=-1, etc.
+    Nk = 2; Nb = 1; Ns = 1
+    rho1 = zeros(ComplexF64, Nk, Nb, Ns, Nk, Nb, Ns)
+    rho1[1, 1, 1, 1, 1, 1] = 1.0
+    rho1[2, 1, 1, 2, 1, 1] = 1.0
+
+    rho2 = zeros(ComplexF64, Nk, Nb, Ns, Nk, Nb, Ns, Nk, Nb, Ns, Nk, Nb, Ns)
+    # ⟨c†_a c†_b c_b c_a⟩ = 1
+    rho2[1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1] = 1.0
+    # ⟨c†_a c†_b c_a c_b⟩ = -1  (antisymmetry under c swap)
+    rho2[1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 1] = -1.0
+    # ⟨c†_b c†_a c_b c_a⟩ = -1  (antisymmetry under c† swap)
+    rho2[2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1] = -1.0
+    # ⟨c†_b c†_a c_a c_b⟩ = 1
+    rho2[2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1] = 1.0
+
+    @testset "normal → density-density, fermionic" begin
+        rho_dd = two_body_normal_to_density_density(rho2, rho1; statistics = :fermi)
+
+        # ⟨c†_a c_b c†_b c_a⟩ = δ_{b,b}⟨c†_a c_a⟩ - ⟨c†_a c†_b c_b c_a⟩ = 1·1 - 1 = 0
+        # but in rho_dd indices are [k1,b1,s1, k3,b3,s3, k2,b2,s2, k4,b4,s4]
+        # i.e. rho_dd[a,b,b,a] = δ_{b,b}⟨c†_a c_a⟩ - rho2[a,b,b,a]
+        @test rho_dd[1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1] ≈ 0.0 atol = 1.0e-14
+
+        # ⟨c†_a c_a c†_b c_b⟩ = δ_{a,a}⟨c†_b... wait, let's use the formula directly:
+        # rho_dd[a,a,b,b] = δ_{a,b}⟨c†_... ⟩ - rho2[a,b,a,b]
+        # a≠b so δ=0: rho_dd[a,a,b,b] = -rho2[a,b,a,b] = -(-1) = 1
+        @test rho_dd[1, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1] ≈ 1.0 atol = 1.0e-14
+    end
+
+    @testset "round-trip: normal → dd → normal" begin
+        rho2_rt = density_density_to_two_body_normal(
+            two_body_normal_to_density_density(rho2, rho1; statistics = :fermi),
+            rho1; statistics = :fermi
+        )
+        @test rho2_rt ≈ rho2 atol = 1.0e-14
+    end
+
+    @testset "round-trip: dd → normal → dd" begin
+        rho_dd = two_body_normal_to_density_density(rho2, rho1; statistics = :fermi)
+        rho_dd_rt = two_body_normal_to_density_density(
+            density_density_to_two_body_normal(rho_dd, rho1; statistics = :fermi),
+            rho1; statistics = :fermi
+        )
+        @test rho_dd_rt ≈ rho_dd atol = 1.0e-14
+    end
+
+    @testset "bosonic sign convention" begin
+        # For bosons: ⟨b†_1 b_3 b†_2 b_4⟩ = δ_{23}⟨b†_1 b_4⟩ + ⟨b†_1 b†_2 b_3 b_4⟩
+        # Use a simple non-antisymmetric rho2 for bosons
+        rho2_bose = zeros(Float64, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1)
+        rho1_bose = zeros(Float64, 2, 1, 1, 2, 1, 1)
+        rho2_bose[1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1] = 2.0
+        rho1_bose[1, 1, 1, 1, 1, 1] = 1.5
+        rho1_bose[2, 1, 1, 2, 1, 1] = 0.5
+
+        rho_dd_b = two_body_normal_to_density_density(rho2_bose, rho1_bose; statistics = :bose)
+        # rho_dd[a,b,b,a] = δ_{b,b}⟨b†_a b_a⟩ + rho2[a,b,b,a] = 1.5 + 2.0 = 3.5
+        @test rho_dd_b[1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1] ≈ 3.5 atol = 1.0e-14
+    end
+end
diff --git a/test/Project.toml b/test/Project.toml
index 8a2b3b0..86e44c3 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -10,6 +10,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 [compat]
 Aqua = "0.8"
 DelimitedFiles = "1.9"
-DiagHamInterface = "1.0,1.1"
+DiagHamInterface = "1.0,1.1,1.2"
 SafeTestsets = "0.1"
 Suppressor = "0.2"
diff --git a/test/test_basics.jl b/test/test_basics.jl
index a8cd196..ac6a309 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -97,7 +97,7 @@ using Test
         # Test normal range numbers
         result = format_with_precision(1.5)
         @test !occursin("e", result) && !occursin("E", result)
-        
+
         # Test negative numbers
         result_neg = format_with_precision(-1.5)
         @test startswith(result_neg, "-")