"""
多粒子系の分配関数（Partition Function）の積の性質を検証するモジュール。
概要:
複数の独立な粒子が存在する系において、全状態和としての分配関数 $Z_{tot}$ が、
各粒子の分配関数 $Z_p$ の積に等しいことを数値的に検証します。
詳細説明:
1. 各粒子に対して、ランダムなエネルギー準位 $\epsilon_{i}$ を生成します。
2. 手法1（全組み合わせ）:
itertools.product を用いて、すべての粒子の状態の組み合わせを網羅し、
全エネルギー $E = \sum \epsilon$ に対するボルツマン因子 $\exp(-E/k_B T)$ を足し合わせます。
3. 手法2（分配関数の積）:
各粒子ごとに分配関数 $Z_p = \sum \exp(-\epsilon/k_B T)$ を計算し、それらを掛け合わせます。
4. 両者の結果が一致することを確認し、独立系における計算の簡略化を理解します。
理論式:
$$Z_{total} = \sum_{s_1, s_2, \dots} e^{-\beta(\epsilon_{s_1} + \epsilon_{s_2} + \dots)} = \prod_{p} \left( \sum_{s} e^{-\beta \epsilon_{s,p}} \right)$$
関連リンク: :doc:`state_sum_usage`
"""
import sys
from numpy import exp
import random
import itertools
[ドキュメント]
def generate_energy_levels(ns):
"""
0から1の範囲でランダムなエネルギー準位（$k_B T$ で正規化済み）を生成します。
"""
return [random.random() for _ in range(ns)]
[ドキュメント]
def main():
"""
分配関数の2通りの計算と、結果の比較を行います。
"""
# ------------------------
# パラメータ設定
# ------------------------
n_particles = 2
n_states = 3
argv = sys.argv
if len(argv) > 1: n_particles = int(argv[1])
if len(argv) > 2: n_states = int(argv[2])
print(f"\n--- Partition Function Verification ---")
print(f"Number of particles (np): {n_particles}")
print(f"Number of states per particle (ns): {n_states}")
# 各粒子のエネルギー準位を生成
eps_list = [generate_energy_levels(n_states) for _ in range(n_particles)]
print("\nGenerated Energy Levels (Normalized by kBT):")
for ip, eps in enumerate(eps_list):
print(f" Particle {ip}: {['{:.4f}'.format(e) for e in eps]}")
# ---------------------------------------------------------
# 手法1: すべての状態の組み合わせを網羅 (Total Combination)
# ---------------------------------------------------------
print("\n[Method 1] Calculating Z_tot using all combinations of states...")
# 状態インデックスの組み合わせを生成
state_indices = [range(n_states)] * n_particles
all_combinations = list(itertools.product(*state_indices))
z_total_comb = 0.0
for ic, combo in enumerate(all_combinations):
# 系の全エネルギーを計算
e_total = sum(eps_list[ip][idx] for ip, idx in enumerate(combo))
# ボルツマン因子の加算
z_contribution = exp(-e_total)
z_total_comb += z_contribution
# 組み合わせが少ない場合のみ詳細を表示
if len(all_combinations) <= 27:
print(f" Combo {ic:2d} {combo}: Etot={e_total:8.5f}, Z_contribution={z_contribution:8.5f}")
print(f"==> Z_total (Method 1) = {z_total_comb:.10f}")
# ---------------------------------------------------------
# 手法2: 各粒子の分配関数の積 (Product of Zp)
# ---------------------------------------------------------
print("\n[Method 2] Calculating Z_tot as a product of individual Zp...")
z_total_prod = 1.0
for ip in range(n_particles):
# 個々の粒子の分配関数 Zp を計算
zp = sum(exp(-e) for e in eps_list[ip])
print(f" Particle {ip} Zp = {zp:.10f}")
z_total_prod *= zp
print(f"==> Z_total (Method 2) = {z_total_prod:.10f}")
# ------------------------
# 結果の比較
# ------------------------
diff = abs(z_total_comb - z_total_prod)
print(f"\nDifference: {diff:.2e}")
if diff < 1e-12:
print("Verification Successful: Z_total = Π Zp")
else:
print("Verification Failed: Check the implementation.")
if __name__ == "__main__":
main()