From 63ca6bf3f0f7db2de36cf88c1d86441b483a6961 Mon Sep 17 00:00:00 2001 From: AFeuerpfeil Date: Sat, 7 Mar 2026 18:52:01 -0500 Subject: [PATCH 1/6] improve printing and bump to 1.1.1 --- Project.toml | 6 +--- src/DiagHamInterface.jl | 2 -- src/utility/numbers.jl | 66 ++++++++++++++++------------------------- 3 files changed, 27 insertions(+), 47 deletions(-) diff --git a/Project.toml b/Project.toml index cffc511..c721a55 100644 --- a/Project.toml +++ b/Project.toml @@ -1,21 +1,17 @@ name = "DiagHamInterface" uuid = "2f27495d-c82b-46d5-b81b-1bb9aa58c4d5" -version = "1.1.0" +version = "1.1.1" authors = ["Andreas Feuerpfeil "] [deps] DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" -Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8" Preferences = "21216c6a-2e73-6563-6e65-726566657250" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -ThreadSafeDicts = "4239201d-c60e-5e0a-9702-85d713665ba7" [compat] DelimitedFiles = "1" -Format = "1.3.7" Preferences = "1.5.0" Printf = "1" SparseArrays = "1" -ThreadSafeDicts = "0.1.6" julia = "1.10" diff --git a/src/DiagHamInterface.jl b/src/DiagHamInterface.jl index 6398097..33ab05f 100644 --- a/src/DiagHamInterface.jl +++ b/src/DiagHamInterface.jl @@ -12,10 +12,8 @@ export write_to_txt export write_matrix_elements using DelimitedFiles -using Format using SparseArrays using Preferences -using ThreadSafeDicts include("utility/backup.jl") include("utility/diagham_path.jl") diff --git a/src/utility/numbers.jl b/src/utility/numbers.jl index 1de0b90..0e3943f 100644 --- a/src/utility/numbers.jl +++ b/src/utility/numbers.jl @@ -1,40 +1,26 @@ -# Maps the format string (e.g. "%.5f") to the compiled formatter Function. -const LOCAL_FORMAT_CACHE = ThreadSafeDict{String, Function}() - """ - format_with_precision(x; atol=eps(float(T)), mode=:auto, maxdigits=20) + format_with_precision(x; atol=0.0, maxdigits=typemax(Int)) -Format number `x` as a string with absolute precision `atol`. -Mode `:auto` uses `%f` for `[1e-3, 1e6)`, else `%e`. Use `:f` or `:e` to force format. +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 = eps(float(T)), mode::Symbol = :auto, maxdigits::Int = 20) where {T <: Real} +function format_with_precision(x::T; atol = 0.0, maxdigits::Int = typemax(Int)) where {T <: Real} iszero(x) && return "0.0" - absx = abs(x) - abs_atol = abs(atol) - - use_e = mode == :e || (mode == :auto && (absx < 1.0e-3 || absx >= 1.0e6)) - - if use_e - exp10 = floor(Int, log10(absx)) - log10_atol = log10(abs_atol) - p = ceil(Int, exp10 - log10_atol) - p = clamp(p, 0, maxdigits) - - formatter = get_threadsafe_formatter("%.$(p)e") - return formatter(x) - else - p = ceil(Int, -log10(abs_atol)) - p = clamp(p, 0, maxdigits) - - formatter = get_threadsafe_formatter("%.$(p)f") - return formatter(x) - end -end - -function get_threadsafe_formatter(fmtstr::String) - return get!(LOCAL_FORMAT_CACHE, fmtstr) do - Format.generate_formatter(fmtstr) + + if iszero(atol) + maxdigits == typemax(Int) && return string(x) + return string(round(x, sigdigits=maxdigits)) end + + mag_x = floor(Int, log10(abs(x))) + mag_atol = floor(Int, log10(abs(atol))) + + required_sigdigits = max(1, mag_x - mag_atol + 1) + s = min(required_sigdigits, maxdigits) + + rounded_x = round(x, sigdigits=s) + return string(rounded_x) end read_number(V::Number) = V @@ -64,18 +50,18 @@ function _read_number_bracket(V::AbstractString) return ComplexF64(real, imag) end -function write_number(V::Number; atol::Real = eps(real(float(V)))) - return format_with_precision(V; atol = atol) +function write_number(V::Number; kwargs...) + return format_with_precision(V; kwargs...) end -function write_number(V::Complex; atol::Real = eps(real(float(V)))) - return string("(", format_with_precision(real(V); atol = atol), ",", format_with_precision(imag(V); atol = atol), ")") +function write_number(V::Complex; kwargs...) + return string("(", format_with_precision(real(V); kwargs...), ",", format_with_precision(imag(V); kwargs...), ")") end -function write_number_space(V::Number; atol::Real = eps(real(float(V)))) - return format_with_precision(V; atol = atol) +function write_number_space(V::Number; kwargs...) + return format_with_precision(V; kwargs...) end -function write_number_space(V::Complex; atol::Real = eps(real(float(V)))) - return "$(format_with_precision(real(V); atol = atol)) $(format_with_precision(imag(V); atol = atol))" +function write_number_space(V::Complex; kwargs...) + return "$(format_with_precision(real(V); kwargs...)) $(format_with_precision(imag(V); kwargs...))" end From 38eaca60313068e334460ce9ababe3a098300015 Mon Sep 17 00:00:00 2001 From: AFeuerpfeil Date: Sat, 7 Mar 2026 18:56:20 -0500 Subject: [PATCH 2/6] fix and format --- src/utility/numbers.jl | 12 ++++++------ test/test_basics.jl | 9 +-------- 2 files changed, 7 insertions(+), 14 deletions(-) diff --git a/src/utility/numbers.jl b/src/utility/numbers.jl index 0e3943f..4907cff 100644 --- a/src/utility/numbers.jl +++ b/src/utility/numbers.jl @@ -7,19 +7,19 @@ 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} iszero(x) && return "0.0" - + if iszero(atol) maxdigits == typemax(Int) && return string(x) - return string(round(x, sigdigits=maxdigits)) + return string(round(x, sigdigits = maxdigits)) end - + mag_x = floor(Int, log10(abs(x))) mag_atol = floor(Int, log10(abs(atol))) - + required_sigdigits = max(1, mag_x - mag_atol + 1) s = min(required_sigdigits, maxdigits) - - rounded_x = round(x, sigdigits=s) + + rounded_x = round(x, sigdigits = s) return string(rounded_x) end diff --git a/test/test_basics.jl b/test/test_basics.jl index 5af2691..a8cd196 100644 --- a/test/test_basics.jl +++ b/test/test_basics.jl @@ -97,14 +97,7 @@ using Test # Test normal range numbers result = format_with_precision(1.5) @test !occursin("e", result) && !occursin("E", result) - - # Test with forced mode - result_e = format_with_precision(1.5; mode = :e) - @test occursin("e", result_e) || occursin("E", result_e) - - result_f = format_with_precision(1.5; mode = :f) - @test !occursin("e", result_f) && !occursin("E", result_f) - + # Test negative numbers result_neg = format_with_precision(-1.5) @test startswith(result_neg, "-") From 0d19acad9ffca3f1f7c49b40b357ac5faa9379a5 Mon Sep 17 00:00:00 2001 From: AFeuerpfeil Date: Thu, 16 Apr 2026 20:55:09 -0400 Subject: [PATCH 3/6] fix multiband ordering of labels in write_matrix_elements --- src/MatrixElements/read.jl | 8 ++++++-- src/MatrixElements/write.jl | 19 +++++++++++++++++++ 2 files changed, 25 insertions(+), 2 deletions(-) diff --git a/src/MatrixElements/read.jl b/src/MatrixElements/read.jl index c0ba567..3493f76 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..0145d62 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") From 3048de05734ab39b40d5d0d214b232c373af07bd Mon Sep 17 00:00:00 2001 From: AFeuerpfeil Date: Fri, 24 Apr 2026 22:16:01 +0200 Subject: [PATCH 4/6] add support for FTIDensity --- Project.toml | 2 +- docs/Project.toml | 2 +- src/DiagHamInterface.jl | 6 + src/FTIDensity/read_density.jl | 263 +++++++++++++++++++++++++++ test/FTIDensity/test_read_density.jl | 247 +++++++++++++++++++++++++ test/Project.toml | 2 +- 6 files changed, 519 insertions(+), 3 deletions(-) create mode 100644 src/FTIDensity/read_density.jl create mode 100644 test/FTIDensity/test_read_density.jl 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 "] [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..c0d6062 --- /dev/null +++ b/src/FTIDensity/read_density.jl @@ -0,0 +1,263 @@ +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 + 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 = 1e-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 = 1e-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/test/FTIDensity/test_read_density.jl b/test/FTIDensity/test_read_density.jl new file mode 100644 index 0000000..348421c --- /dev/null +++ b/test/FTIDensity/test_read_density.jl @@ -0,0 +1,247 @@ +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 \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' \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' \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 \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 +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' +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' +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' +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 \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 \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 +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=1e-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=1e-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=1e-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=1e-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=1e-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" From cc3602e3ff0e43a9e994e714ecc9f6559101bc4d Mon Sep 17 00:00:00 2001 From: AFeuerpfeil Date: Fri, 24 Apr 2026 22:17:47 +0200 Subject: [PATCH 5/6] restore default --- src/utility/numbers.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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) From ce0f419f872f8d19b0b2d64948d3739410c97308 Mon Sep 17 00:00:00 2001 From: AFeuerpfeil Date: Fri, 24 Apr 2026 22:19:51 +0200 Subject: [PATCH 6/6] format --- src/FTIDensity/read_density.jl | 34 +++++++------ src/MatrixElements/read.jl | 2 +- src/MatrixElements/write.jl | 2 +- test/FTIDensity/test_read_density.jl | 76 ++++++++++++++-------------- test/test_basics.jl | 2 +- 5 files changed, 60 insertions(+), 56 deletions(-) diff --git a/src/FTIDensity/read_density.jl b/src/FTIDensity/read_density.jl index c0d6062..f91956f 100644 --- a/src/FTIDensity/read_density.jl +++ b/src/FTIDensity/read_density.jl @@ -59,7 +59,7 @@ end function _determine_type(values::AbstractVector, atol::Real) isempty(values) && return Float64 - any(v -> abs(imag(v)) > atol, values) ? ComplexF64 : Float64 + return any(v -> abs(imag(v)) > atol, values) ? ComplexF64 : Float64 end """ @@ -78,7 +78,7 @@ data. For single-band files Nband=1 and Nspin=1. For files with a `spin` column Since FTIDensity only writes diagonal-in-k entries, all off-diagonal-k entries are zero. """ -function read_one_body_density(filename::AbstractString; atol::Real = 1e-14) +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") @@ -95,7 +95,7 @@ function read_one_body_density(filename::AbstractString; atol::Real = 1e-14) ky_col = index_matrix[:, 2] Nkx = maximum(kx_col) + 1 Nky = maximum(ky_col) + 1 - Nk = Nkx * Nky + Nk = Nkx * Nky if Ncols_index == 2 # 1-band: kx ky val @@ -115,7 +115,7 @@ function read_one_body_density(filename::AbstractString; atol::Real = 1e-14) 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 + 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) @@ -123,15 +123,15 @@ function read_one_body_density(filename::AbstractString; atol::Real = 1e-14) end else # 4/6-band: kx ky spin sigma_cr sigma_an val - spin_col = index_matrix[:, 3] + 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 + 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) @@ -152,7 +152,7 @@ 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 = 1e-14) +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") @@ -168,17 +168,17 @@ function read_two_body_density(filename::AbstractString; atol::Real = 1e-14) 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] + 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] + 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 + 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 @@ -188,7 +188,7 @@ function read_two_body_density(filename::AbstractString; atol::Real = 1e-14) 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) + 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 @@ -213,7 +213,8 @@ For bosons (`statistics=:bose`, commutator `[b,b†]=δ`): 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} + 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]) @@ -235,7 +236,8 @@ 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} + 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]) diff --git a/src/MatrixElements/read.jl b/src/MatrixElements/read.jl index 3493f76..14252a2 100644 --- a/src/MatrixElements/read.jl +++ b/src/MatrixElements/read.jl @@ -60,7 +60,7 @@ function reinstate_indices(header, indices, coeffs) ## This function takes care label_priority = h -> begin base = split(h, "_")[1] startswith(base, standard_band_label()) ? (0, base) : - startswith(base, standard_valley_label()) ? (1, base) : (2, base) + startswith(base, standard_valley_label()) ? (1, base) : (2, base) end local_order = sortperm(header[pos], by = label_priority) append!(perm, pos[local_order]) diff --git a/src/MatrixElements/write.jl b/src/MatrixElements/write.jl index 0145d62..2eab881 100644 --- a/src/MatrixElements/write.jl +++ b/src/MatrixElements/write.jl @@ -78,7 +78,7 @@ 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) + startswith(base, standard_valley_label()) ? (1, base) : (2, base) end suffixes = [parse(Int, split(h, "_")[2]) for h in label] perm = Int[] diff --git a/test/FTIDensity/test_read_density.jl b/test/FTIDensity/test_read_density.jl index 348421c..8b14f9a 100644 --- a/test/FTIDensity/test_read_density.jl +++ b/test/FTIDensity/test_read_density.jl @@ -52,13 +52,13 @@ end @testset "read_one_body_density" begin @testset "1-band, 2x2 lattice, real" begin content = """# kx ky -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 -""" + 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) @@ -78,11 +78,11 @@ end @testset "2-band, 1x1 lattice, complex" begin content = """# kx ky sigma sigma' -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 -""" + 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) @@ -100,15 +100,15 @@ end @testset "4-band (spin), 1x1 lattice" begin content = """# kx ky spin sigma sigma' -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 -""" + 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) @@ -128,12 +128,12 @@ end @testset "Comment lines interspersed with data" begin # DiagHam inserts partial-density comments between data blocks content = """# kx ky sigma sigma' -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) -""" + 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) @@ -163,9 +163,9 @@ end # 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 -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 -""" + 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) @@ -208,27 +208,29 @@ end # ⟨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=1e-14 + @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=1e-14 + @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=1e-14 + 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=1e-14 + rho1; statistics = :fermi + ) + @test rho_dd_rt ≈ rho_dd atol = 1.0e-14 end @testset "bosonic sign convention" begin @@ -242,6 +244,6 @@ end 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=1e-14 + @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/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, "-")