25113611074

Committed 29 Apr 2026 02:04PM UTC coverage: 84.91% (+1.0%) from 83.942%

Build # 25113611074

Build Type

Pull #129

github

Committed by

web-flow

Commit Message

Merge ae5542e01 into 1a743d79c

Pull Request Pull Request #129: Grand-canonical nested sampling for lattice systems

Coverage Stats

207 of 215 new or added lines in 3 files covered. (96.28%)

3 existing lines in 1 file now uncovered.

2065 of 2432 relevant lines covered (84.91%)

615272.14 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

96.67

/src/AnalysisTools/AnalysisTools.jl

"""
    AnalysisTools

Module for analyzing the output of the sampling.
"""
module AnalysisTools

using DataFrames
using CSV, Arrow

export read_output
export ωᵢ, partition_function, internal_energy, cv
export gc_thermodynamic_stats

"""
    read_output(filename::String)

Reads the output file and returns a DataFrame.
"""
function read_output(filename::String)
    if splitext(filename)[end] == ".csv"
        data = CSV.File(filename)
    elseif splitext(filename)[end] == ".arrow"
        data = Arrow.Table(filename)
    else
        error("Unsupported file format. Please provide a .csv or .arrow file.")
    end
    return DataFrame(data)
end

"""
    ωᵢ(iters::Vector{Int}, n_walkers::Int; n_cull::Int=1, ω0::Float64=1.0)

Calculates the \$\\omega\$ factors for the given number of iterations and walkers.
The \$\\omega\$ factors account for the fractions of phase-space volume sampled during
each nested sampling iteration, defined as:
```math
\\omega_i = \\frac{C}{K+C} \\left(\\frac{K}{K+C}\\right)^i
```
where \$K\$ is the number of walkers, \$C\$ is the number of culled walkers, 
and \$i\$ is the iteration number.

# Arguments
- `iters::Vector{Int}`: The iteration numbers.
- `n_walkers::Int`: The number of walkers.
- `n_cull::Int`: The number of culled walkers. Default is 1.
- `ω0::Float64`: The initial \$\\omega\$ factor. Default is 1.0.

# Returns
- A vector of \$\\omega\$ factors.
"""
function ωᵢ(iters::AbstractVector{Int}, n_walkers::Int; n_cull::Int=1, ω0::Float64=1.0)
    ωi = ω0 * (n_cull/(n_walkers+n_cull)) * (n_walkers/(n_walkers+n_cull)).^iters
    return ωi
end

"""
    partition_function(β::Float64, ωi::Vector{Float64}, Ei::Vector{Float64})

Calculates the partition function for the given \$\\beta\$, \$\\omega\$ factors, and energies.
The partition function is defined as:
```math
Z(\\beta) = \\sum_i \\omega_i \\exp(-E_i \\beta)
```
where \$\\omega_i\$ is the \$i\$-th \$\\omega\$ factor, \$E_i\$ is the \$i\$-th energy, and \$\\beta\$ is the inverse temperature.

# Arguments
- `β::Float64`: The inverse temperature.
- `ωi::Vector{Float64}`: The \$\\omega\$ factors.
- `Ei::Vector{Float64}`: The energies.

# Returns
- The partition function.
"""
function partition_function(β::Float64, 
                            ωi::Vector{Float64}, 
                            Ei::Vector{Float64})
    z = sum(ωi.*exp.(-Ei.*β))
    return z
end

"""
    internal_energy(β::Float64, ωi::Vector{Float64}, ei::Vector{Float64})

Calculates the internal energy from the partition function for the given \$\\beta\$, \$\\omega\$ factors, and energies.
The internal energy is defined as:
```math
U(\\beta) = \\frac{\\sum_i \\omega_i E_i \\exp(-E_i \\beta)}{\\sum_i \\omega_i \\exp(-E_i \\beta)}
```
where \$\\omega_i\$ is the \$i\$-th \$\\omega\$ factor, \$E_i\$ is the \$i\$-th energy, and \$\\beta\$ is the inverse temperature.

# Arguments
- `β::Float64`: The inverse temperature.
- `ωi::Vector{Float64}`: The \$\\omega\$ factors.
- `Ei::Vector{Float64}`: The energies in eV.

# Returns
- The internal energy.
"""
function internal_energy(β::Float64, 
                         ωi::Vector{Float64}, 
                         Ei::Vector{Float64})
    u = sum(ωi.*Ei.*exp.(-Ei.*β))/sum(ωi.*exp.(-Ei.*β))
    return u
end

"""
    cv(β::Float64, omega_i::Vector{Float64}, Ei::Vector{Float64}, dof::Int)

Calculates the constant-volume heat capacity for the given \$\\beta\$, \$\\omega\$ factors, energies, and degrees of freedom.
The heat capacity is defined as:
```math
C_V(\\beta) = \\frac{\\mathrm{dof} \\cdot k_B}{2} + k_B \\beta^2 \\left(\\frac{\\sum_i \\omega_i E_i^2 \\exp(-E_i \\beta)}{Z(\\beta)} - U(\\beta)^2\\right)
```
where \$\\mathrm{dof}\$ is the degrees of freedom, \$k_B\$ is the Boltzmann constant (in units of eV/K), \$\\beta\$ is the inverse temperature, 
\$\\omega_i\$ is the \$i\$-th \$\\omega\$ factor, \$E_i\$ is the \$i\$-th energy, \$Z(\\beta)\$ is the partition function, and \$U(\\beta)\$ is the internal energy.

# Arguments
- `β::Float64`: The inverse temperature.
- `ωi::Vector{Float64}`: The \$\\omega\$ factors.
- `Ei::Vector{Float64}`: The energies in eV.
- `dof::Int`: The degrees of freedom, equals to the number of dimensions times the number of particles.

# Returns
- The constant-volume heat capacity.
"""
function cv(β::Float64,
            ωi::Vector{Float64}, 
            Ei::Vector{Float64},
            dof::Int64;
            kb::Float64=8.617333262e-5)
    expo = ωi.*exp.(-Ei.*β)
    ei_expo = Ei.*expo
    ei2_expo = Ei.*ei_expo
    z = sum(expo)
    u = sum(ei_expo)/z
    cv = dof*kb/2.0 + kb*β^2 * (sum(ei2_expo)/z - u^2)
    return cv
end

"""
    cv(df::DataFrame, βs::Vector{Float64}, dof::Int, n_walkers::Int)

(Nested Sampling) Calculates the constant-volume heat capacity at constant volume for the given DataFrame, inverse temperatures, degrees of freedom, and number of walkers.
The heat capacity is defined as:
```math
C_V(\\beta) = \\frac{\\mathrm{dof} \\cdot k_B}{2} + k_B \\beta^2 \\left(\\frac{\\sum_i \\omega_i E_i^2 \\exp(-E_i \\beta)}{Z(\\beta)} - U(\\beta)^2\\right)
```
where \$\\mathrm{dof}\$ is the degrees of freedom, \$k_B\$ is the Boltzmann constant (in units of eV/K), \$\\beta\$ is the inverse temperature,
\$\\omega_i\$ is the \$i\$-th \$\\omega\$ factor, \$E_i\$ is the \$i\$-th energy, \$Z(\\beta)\$ is the partition function, and \$U(\\beta)\$ is the internal energy.

# Arguments
- `df::DataFrame`: The DataFrame containing the output data.
- `βs::Vector{Float64}`: Inverse temperatures.
- `dof::Int`: The degrees of freedom, equal to the number of dimensions times the number of particles. For a lattice, it is zero.
- `n_walkers::Int`: The number of walkers.
- `n_cull::Int`: The number of culled walkers. Default is 1.
- `ω0::Float64`: The initial \$\\omega\$ factor. Default is 1.0.

# Returns
- A vector of constant-volume heat capacities.
"""
function cv(df::DataFrame, 
            βs::Vector{Float64}, 
            dof::Int, n_walkers::Int; 
            n_cull::Int=1, 
            ω0::Float64=1.0, 
            kb::Float64=8.617333262e-5)
    ωi = ωᵢ(df.iter, n_walkers; n_cull=n_cull, ω0=ω0)
    Ei = df.emax .- minimum(df.emax)
    cvs = Vector{Float64}(undef, length(βs))
    Threads.@threads for (i, b) in collect(enumerate(βs))
        cvs[i] = cv(b, ωi, Ei, dof; kb=kb)
    end
    return cvs
end

"""
    cv(Ts::Vector{Float64}, dof::Int, energy_bins::Vector{Float64}, entropy::Vector{Float64})

(Wang-Landau Sampling) Calculates the constant-volume heat capacity at constant volume for the given temperatures, degrees of 
freedom, energy bins, and entropy. The kinetic energy is treated classically, and is added to the heat capacity as \$dof \\cdot k_B/2\$.

# Arguments
- `Ts::Vector{Float64}`: The temperatures in Kelvin.
- `dof::Int`: The degrees of freedom, equals to the number of dimensions times the number of particles. For a lattice, it is zero.
- `energy_bins::Vector{Float64}`: The energy bins in eV.
- `entropy::Vector{Float64}`: The entropy.

# Returns
- A vector of constant-volume heat capacities.
"""
function cv(Ts::Vector{Float64}, dof::Int, energy_bins::Vector{Float64}, entropy::Vector{Float64})
    kb = 8.617333262e-5 # eV/K
    β = 1 ./(kb.*Ts)
    E_rel = energy_bins .- minimum(energy_bins)
    S_shifted = entropy .- minimum(entropy[entropy .> 0])
    g = exp.(S_shifted)
    Z = zeros(length(Ts))
    E_avg = zeros(length(Ts))
    E2_avg = zeros(length(Ts))
    Cv = zeros(length(Ts))
    for (i, temp) in enumerate(Ts)
        Z[i] = sum(exp.(-E_rel ./ (kb * temp)) .* g)
        E_avg[i] = sum(E_rel .* exp.(-E_rel ./ (kb * temp)) .* g) / Z[i]
        E2_avg[i] = sum(E_rel.^2 .* exp.(-E_rel ./ (kb * temp)) .* g) / Z[i]
        Cv[i] = (E2_avg[i] .- E_avg[i].^2) ./ (kb * temp.^2) .+ dof*kb/2
    end
    return Cv
end


"""
    gc_thermodynamic_stats(β::Float64, ωi::Vector{Float64},
                           grand_energies::Vector{Float64},
                           energies::Vector{Float64},
                           numbers::Vector{Int},
                           μ::Float64;
                           kb::Float64=8.617333262e-5)

Compute grand-canonical thermodynamic averages from nested sampling output.

The log-sum-exp trick is used for numerical stability. The grand-canonical
heat capacity at constant μ is:

    C_{V,μ} = k_B β² [Var(E) − μ Cov(E, N)]

# Arguments
- `β::Float64`: Inverse temperature 1/(k_B T).
- `ωi::Vector{Float64}`: Phase-space volume weights from NS.
- `grand_energies::Vector{Float64}`: Ω_i = E_i − μ N_i values.
- `energies::Vector{Float64}`: E_i values.
- `numbers::Vector{Int}`: N_i values.
- `μ::Float64`: Chemical potential.
- `kb::Float64`: Boltzmann constant (default: eV/K).

# Returns
- `(⟨E⟩, C_{V,μ}, ⟨N⟩)`: Mean energy, GC heat capacity, mean particle number.
"""
function gc_thermodynamic_stats(β::Float64,
                                 ωi::Vector{Float64},
                                 grand_energies::Vector{Float64},
                                 energies::Vector{Float64},
                                 numbers::Vector{Int},
                                 μ::Float64;
                                 kb::Float64=8.617333262e-5)
    n = length(ωi)
    if n != length(grand_energies) || n != length(energies) || n != length(numbers)
        throw(DimensionMismatch("All input vectors must have the same length"))
    end
    if n == 0
        return NaN, NaN, NaN
    end

    # Log-sum-exp for numerical stability
    log_terms = [log(ωi[i]) - β * grand_energies[i] for i in 1:n]
    max_log = maximum(log_terms)

    z = 0.0
    u = 0.0   # ⟨E⟩
    u2 = 0.0  # ⟨E²⟩
    n_sum = 0.0  # ⟨N⟩
    en_sum = 0.0 # ⟨EN⟩

    for i in 1:n
        w = exp(log_terms[i] - max_log)
        z += w
        u += w * energies[i]
        u2 += w * energies[i]^2
        n_sum += w * numbers[i]
        en_sum += w * energies[i] * numbers[i]
    end

    if z == 0.0
        return NaN, NaN, NaN
    end

    u /= z
    u2 /= z
    n_avg = n_sum / z
    en_avg = en_sum / z

    var_e = u2 - u^2
    cov_en = en_avg - u * n_avg

    cv = kb * β^2 * (var_e - μ * cov_en)

    return u, cv, n_avg
end

"""
    gc_thermodynamic_stats(df::DataFrame, βs::Vector{Float64},
                           n_walkers::Int, μ::Float64;
                           n_cull::Int=1, ω0::Float64=1.0,
                           kb::Float64=8.617333262e-5)

Compute grand-canonical thermodynamic stats from a GC-NS output DataFrame.

The DataFrame must have columns `:iter`, `:omega`, `:energy`, `:num_particles`.
Adds the live-walker contribution to the end of the recorded samples for correct
normalization.

# Arguments
- `df::DataFrame`: GC-NS output with columns `[:iter, :omega, :energy, :num_particles]`.
- `βs::Vector{Float64}`: Inverse temperatures at which to evaluate.
- `n_walkers::Int`: Number of walkers used in the NS run.
- `μ::Float64`: Chemical potential.
- `n_cull::Int=1`: Number of walkers culled per iteration.
- `ω0::Float64=1.0`: Initial phase-space volume.
- `kb::Float64`: Boltzmann constant (default: eV/K).

# Returns
- `(mean_E, Cv, mean_N)`: Vectors of ⟨E⟩, C_{V,μ}, and ⟨N⟩ at each β.
"""
function gc_thermodynamic_stats(df::DataFrame,
                                 βs::Vector{Float64},
                                 n_walkers::Int,
                                 μ::Float64;
                                 n_cull::Int=1,
                                 ω0::Float64=1.0,
                                 kb::Float64=8.617333262e-5)
    ωi = ωᵢ(df.iter, n_walkers; n_cull=n_cull, ω0=ω0)
    grand_es = df.omega
    Es = df.energy
    Ns = df.num_particles

    mean_Es = Vector{Float64}(undef, length(βs))
    Cvs = Vector{Float64}(undef, length(βs))
    mean_Ns = Vector{Float64}(undef, length(βs))

    Threads.@threads for (i, b) in collect(enumerate(βs))
        mean_Es[i], Cvs[i], mean_Ns[i] = gc_thermodynamic_stats(
            b, ωi, grand_es, Es, Ns, μ; kb=kb)
    end

    return mean_Es, Cvs, mean_Ns
end


end # module AnalysisTools

1	"""
2	AnalysisTools
3
4	Module for analyzing the output of the sampling.
5	"""
6	module AnalysisTools
7
8	using DataFrames
9	using CSV, Arrow
10
11	export read_output
12	export ωᵢ, partition_function, internal_energy, cv
13	export gc_thermodynamic_stats
14
15	"""
16	read_output(filename::String)
17
18	Reads the output file and returns a DataFrame.
19	"""
20	function read_output(filename::String)	8✔
21	if splitext(filename)[end] == ".csv"	8✔
22	data = CSV.File(filename)	4✔
23	elseif splitext(filename)[end] == ".arrow"	4✔
24	data = Arrow.Table(filename)	4✔
25	else
26	error("Unsupported file format. Please provide a .csv or .arrow file.")	×
27	end
28	return DataFrame(data)	8✔
29	end
30
31	"""
32	ωᵢ(iters::Vector{Int}, n_walkers::Int; n_cull::Int=1, ω0::Float64=1.0)
33
34	Calculates the \$\\omega\$ factors for the given number of iterations and walkers.
35	The \$\\omega\$ factors account for the fractions of phase-space volume sampled during
36	each nested sampling iteration, defined as:
37	```math
38	\\omega_i = \\frac{C}{K+C} \\left(\\frac{K}{K+C}\\right)^i
39	```
40	where \$K\$ is the number of walkers, \$C\$ is the number of culled walkers,
41	and \$i\$ is the iteration number.
42
43	# Arguments
44	- `iters::Vector{Int}`: The iteration numbers.
45	- `n_walkers::Int`: The number of walkers.
46	- `n_cull::Int`: The number of culled walkers. Default is 1.
47	- `ω0::Float64`: The initial \$\\omega\$ factor. Default is 1.0.
48
49	# Returns
50	- A vector of \$\\omega\$ factors.
51	"""
52	function ωᵢ(iters::AbstractVector{Int}, n_walkers::Int; n_cull::Int=1, ω0::Float64=1.0)	72✔
53	ωi = ω0 * (n_cull/(n_walkers+n_cull)) * (n_walkers/(n_walkers+n_cull)).^iters	36✔
54	return ωi	36✔
55	end
56
57	"""
58	partition_function(β::Float64, ωi::Vector{Float64}, Ei::Vector{Float64})
59
60	Calculates the partition function for the given \$\\beta\$, \$\\omega\$ factors, and energies.
61	The partition function is defined as:
62	```math
63	Z(\\beta) = \\sum_i \\omega_i \\exp(-E_i \\beta)
64	```
65	where \$\\omega_i\$ is the \$i\$-th \$\\omega\$ factor, \$E_i\$ is the \$i\$-th energy, and \$\\beta\$ is the inverse temperature.
66
67	# Arguments
68	- `β::Float64`: The inverse temperature.
69	- `ωi::Vector{Float64}`: The \$\\omega\$ factors.
70	- `Ei::Vector{Float64}`: The energies.
71
72	# Returns
73	- The partition function.
74	"""
75	function partition_function(β::Float64,	16✔
76	ωi::Vector{Float64},
77	Ei::Vector{Float64})
78	z = sum(ωi.exp.(-Ei.β))	24✔
79	return z	12✔
80	end
81
82	"""
83	internal_energy(β::Float64, ωi::Vector{Float64}, ei::Vector{Float64})
84
85	Calculates the internal energy from the partition function for the given \$\\beta\$, \$\\omega\$ factors, and energies.
86	The internal energy is defined as:
87	```math
88	U(\\beta) = \\frac{\\sum_i \\omega_i E_i \\exp(-E_i \\beta)}{\\sum_i \\omega_i \\exp(-E_i \\beta)}
89	```
90	where \$\\omega_i\$ is the \$i\$-th \$\\omega\$ factor, \$E_i\$ is the \$i\$-th energy, and \$\\beta\$ is the inverse temperature.
91
92	# Arguments
93	- `β::Float64`: The inverse temperature.
94	- `ωi::Vector{Float64}`: The \$\\omega\$ factors.
95	- `Ei::Vector{Float64}`: The energies in eV.
96
97	# Returns
98	- The internal energy.
99	"""
100	function internal_energy(β::Float64,	12✔
101	ωi::Vector{Float64},
102	Ei::Vector{Float64})
103	u = sum(ωi.Ei.exp.(-Ei.β))/sum(ωi.exp.(-Ei.*β))	20✔
104	return u	12✔
105	end
106
107	"""
108	cv(β::Float64, omega_i::Vector{Float64}, Ei::Vector{Float64}, dof::Int)
109
110	Calculates the constant-volume heat capacity for the given \$\\beta\$, \$\\omega\$ factors, energies, and degrees of freedom.
111	The heat capacity is defined as:
112	```math
113	C_V(\\beta) = \\frac{\\mathrm{dof} \\cdot k_B}{2} + k_B \\beta^2 \\left(\\frac{\\sum_i \\omega_i E_i^2 \\exp(-E_i \\beta)}{Z(\\beta)} - U(\\beta)^2\\right)
114	```
115	where \$\\mathrm{dof}\$ is the degrees of freedom, \$k_B\$ is the Boltzmann constant (in units of eV/K), \$\\beta\$ is the inverse temperature,
116	\$\\omega_i\$ is the \$i\$-th \$\\omega\$ factor, \$E_i\$ is the \$i\$-th energy, \$Z(\\beta)\$ is the partition function, and \$U(\\beta)\$ is the internal energy.
117
118	# Arguments
119	- `β::Float64`: The inverse temperature.
120	- `ωi::Vector{Float64}`: The \$\\omega\$ factors.
121	- `Ei::Vector{Float64}`: The energies in eV.
122	- `dof::Int`: The degrees of freedom, equals to the number of dimensions times the number of particles.
123
124	# Returns
125	- The constant-volume heat capacity.
126	"""
127	function cv(β::Float64,	64✔
128	ωi::Vector{Float64},
129	Ei::Vector{Float64},
130	dof::Int64;
131	kb::Float64=8.617333262e-5)
132	expo = ωi.exp.(-Ei.β)	60✔
133	ei_expo = Ei.*expo	56✔
134	ei2_expo = Ei.*ei_expo	56✔
135	z = sum(expo)	52✔
136	u = sum(ei_expo)/z	52✔
137	cv = dofkb/2.0 + kbβ^2 * (sum(ei2_expo)/z - u^2)	52✔
138	return cv	28✔
139	end
140
141	"""
142	cv(df::DataFrame, βs::Vector{Float64}, dof::Int, n_walkers::Int)
143
144	(Nested Sampling) Calculates the constant-volume heat capacity at constant volume for the given DataFrame, inverse temperatures, degrees of freedom, and number of walkers.
145	The heat capacity is defined as:
146	```math
147	C_V(\\beta) = \\frac{\\mathrm{dof} \\cdot k_B}{2} + k_B \\beta^2 \\left(\\frac{\\sum_i \\omega_i E_i^2 \\exp(-E_i \\beta)}{Z(\\beta)} - U(\\beta)^2\\right)
148	```
149	where \$\\mathrm{dof}\$ is the degrees of freedom, \$k_B\$ is the Boltzmann constant (in units of eV/K), \$\\beta\$ is the inverse temperature,
150	\$\\omega_i\$ is the \$i\$-th \$\\omega\$ factor, \$E_i\$ is the \$i\$-th energy, \$Z(\\beta)\$ is the partition function, and \$U(\\beta)\$ is the internal energy.
151
152	# Arguments
153	- `df::DataFrame`: The DataFrame containing the output data.
154	- `βs::Vector{Float64}`: Inverse temperatures.
155	- `dof::Int`: The degrees of freedom, equal to the number of dimensions times the number of particles. For a lattice, it is zero.
156	- `n_walkers::Int`: The number of walkers.
157	- `n_cull::Int`: The number of culled walkers. Default is 1.
158	- `ω0::Float64`: The initial \$\\omega\$ factor. Default is 1.0.
159
160	# Returns
161	- A vector of constant-volume heat capacities.
162	"""
163	function cv(df::DataFrame,	8✔
164	βs::Vector{Float64},
165	dof::Int, n_walkers::Int;
166	n_cull::Int=1,
167	ω0::Float64=1.0,
168	kb::Float64=8.617333262e-5)
169	ωi = ωᵢ(df.iter, n_walkers; n_cull=n_cull, ω0=ω0)	4✔
170	Ei = df.emax .- minimum(df.emax)	4✔
171	cvs = Vector{Float64}(undef, length(βs))	4✔
172	Threads.@threads for (i, b) in collect(enumerate(βs))	8✔
173	cvs[i] = cv(b, ωi, Ei, dof; kb=kb)	12✔
174	end
175	return cvs	4✔
176	end
177
178	"""
179	cv(Ts::Vector{Float64}, dof::Int, energy_bins::Vector{Float64}, entropy::Vector{Float64})
180
181	(Wang-Landau Sampling) Calculates the constant-volume heat capacity at constant volume for the given temperatures, degrees of
182	freedom, energy bins, and entropy. The kinetic energy is treated classically, and is added to the heat capacity as \$dof \\cdot k_B/2\$.
183
184	# Arguments
185	- `Ts::Vector{Float64}`: The temperatures in Kelvin.
186	- `dof::Int`: The degrees of freedom, equals to the number of dimensions times the number of particles. For a lattice, it is zero.
187	- `energy_bins::Vector{Float64}`: The energy bins in eV.
188	- `entropy::Vector{Float64}`: The entropy.
189
190	# Returns
191	- A vector of constant-volume heat capacities.
192	"""
193	function cv(Ts::Vector{Float64}, dof::Int, energy_bins::Vector{Float64}, entropy::Vector{Float64})	4✔
194	kb = 8.617333262e-5 # eV/K	4✔
195	β = 1 ./(kb.*Ts)	8✔
196	E_rel = energy_bins .- minimum(energy_bins)	8✔
197	S_shifted = entropy .- minimum(entropy[entropy .> 0])	8✔
198	g = exp.(S_shifted)	8✔
199	Z = zeros(length(Ts))	8✔
200	E_avg = zeros(length(Ts))	8✔
201	E2_avg = zeros(length(Ts))	8✔
202	Cv = zeros(length(Ts))	8✔
203	for (i, temp) in enumerate(Ts)	4✔
204	Z[i] = sum(exp.(-E_rel ./ (kb * temp)) .* g)	10✔
205	E_avg[i] = sum(E_rel .* exp.(-E_rel ./ (kb * temp)) .* g) / Z[i]	10✔
206	E2_avg[i] = sum(E_rel.^2 .* exp.(-E_rel ./ (kb * temp)) .* g) / Z[i]	10✔
207	Cv[i] = (E2_avg[i] .- E_avg[i].^2) ./ (kb * temp.^2) .+ dof*kb/2	8✔
208	end	12✔
209	return Cv	4✔
210	end
211
212
213	"""
214	gc_thermodynamic_stats(β::Float64, ωi::Vector{Float64},
215	grand_energies::Vector{Float64},
216	energies::Vector{Float64},
217	numbers::Vector{Int},
218	μ::Float64;
219	kb::Float64=8.617333262e-5)
220
221	Compute grand-canonical thermodynamic averages from nested sampling output.
222
223	The log-sum-exp trick is used for numerical stability. The grand-canonical
224	heat capacity at constant μ is:
225
226	C_{V,μ} = k_B β² [Var(E) − μ Cov(E, N)]
227
228	# Arguments
229	- `β::Float64`: Inverse temperature 1/(k_B T).
230	- `ωi::Vector{Float64}`: Phase-space volume weights from NS.
231	- `grand_energies::Vector{Float64}`: Ω_i = E_i − μ N_i values.
232	- `energies::Vector{Float64}`: E_i values.
233	- `numbers::Vector{Int}`: N_i values.
234	- `μ::Float64`: Chemical potential.
235	- `kb::Float64`: Boltzmann constant (default: eV/K).
236
237	# Returns
238	- `(⟨E⟩, C_{V,μ}, ⟨N⟩)`: Mean energy, GC heat capacity, mean particle number.
239	"""
240	function gc_thermodynamic_stats(β::Float64,	32✔
241	ωi::Vector{Float64},
242	grand_energies::Vector{Float64},
243	energies::Vector{Float64},
244	numbers::Vector{Int},
245	μ::Float64;
246	kb::Float64=8.617333262e-5)
247	n = length(ωi)	16✔
248	if n != length(grand_energies) \|\| n != length(energies) \|\| n != length(numbers)	28✔
249	throw(DimensionMismatch("All input vectors must have the same length"))	4✔
250	end
251	if n == 0	12✔
NEW 252	return NaN, NaN, NaN	×
253	end
254
255	# Log-sum-exp for numerical stability
256	log_terms = [log(ωi[i]) - β * grand_energies[i] for i in 1:n]	12✔
257	max_log = maximum(log_terms)	12✔
258
259	z = 0.0	12✔
260	u = 0.0 # ⟨E⟩	12✔
261	u2 = 0.0 # ⟨E²⟩	12✔
262	n_sum = 0.0 # ⟨N⟩	12✔
263	en_sum = 0.0 # ⟨EN⟩	12✔
264
265	for i in 1:n	12✔
266	w = exp(log_terms[i] - max_log)	14,164✔
267	z += w	14,164✔
268	u += w * energies[i]	14,164✔
269	u2 += w * energies[i]^2	14,164✔
270	n_sum += w * numbers[i]	14,164✔
271	en_sum += w * energies[i] * numbers[i]	14,164✔
272	end	28,316✔
273
274	if z == 0.0	12✔
NEW 275	return NaN, NaN, NaN	×
276	end
277
278	u /= z	12✔
279	u2 /= z	12✔
280	n_avg = n_sum / z	12✔
281	en_avg = en_sum / z	12✔
282
283	var_e = u2 - u^2	12✔
284	cov_en = en_avg - u * n_avg	12✔
285
286	cv = kb * β^2 * (var_e - μ * cov_en)	12✔
287
288	return u, cv, n_avg	12✔
289	end
290
291	"""
292	gc_thermodynamic_stats(df::DataFrame, βs::Vector{Float64},
293	n_walkers::Int, μ::Float64;
294	n_cull::Int=1, ω0::Float64=1.0,
295	kb::Float64=8.617333262e-5)
296
297	Compute grand-canonical thermodynamic stats from a GC-NS output DataFrame.
298
299	The DataFrame must have columns `:iter`, `:omega`, `:energy`, `:num_particles`.
300	Adds the live-walker contribution to the end of the recorded samples for correct
301	normalization.
302
303	# Arguments
304	- `df::DataFrame`: GC-NS output with columns `[:iter, :omega, :energy, :num_particles]`.
305	- `βs::Vector{Float64}`: Inverse temperatures at which to evaluate.
306	- `n_walkers::Int`: Number of walkers used in the NS run.
307	- `μ::Float64`: Chemical potential.
308	- `n_cull::Int=1`: Number of walkers culled per iteration.
309	- `ω0::Float64=1.0`: Initial phase-space volume.
310	- `kb::Float64`: Boltzmann constant (default: eV/K).
311
312	# Returns
313	- `(mean_E, Cv, mean_N)`: Vectors of ⟨E⟩, C_{V,μ}, and ⟨N⟩ at each β.
314	"""
315	function gc_thermodynamic_stats(df::DataFrame,	16✔
316	βs::Vector{Float64},
317	n_walkers::Int,
318	μ::Float64;
319	n_cull::Int=1,
320	ω0::Float64=1.0,
321	kb::Float64=8.617333262e-5)
322	ωi = ωᵢ(df.iter, n_walkers; n_cull=n_cull, ω0=ω0)	8✔
323	grand_es = df.omega	8✔
324	Es = df.energy	8✔
325	Ns = df.num_particles	8✔
326
327	mean_Es = Vector{Float64}(undef, length(βs))	8✔
328	Cvs = Vector{Float64}(undef, length(βs))	8✔
329	mean_Ns = Vector{Float64}(undef, length(βs))	8✔
330
331	Threads.@threads for (i, b) in collect(enumerate(βs))	16✔
332	mean_Es[i], Cvs[i], mean_Ns[i] = gc_thermodynamic_stats(	8✔
333	b, ωi, grand_es, Es, Ns, μ; kb=kb)
334	end
335
336	return mean_Es, Cvs, mean_Ns	8✔
337	end
338
339
340	end # module AnalysisTools

wexlergroup / FreeBird.jl / 25113611074

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous