#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Seebeck2.py と Seebeck_ZT.py の計算結果比較テスト。 比較対象: - キャリア濃度 ne - Seebeck 係数 S - Lorentz 数 L - 移動度 mu 比較条件: - r: 既定値 [-0.5, 0, 0.5, 1.0, 1.5, 2.0] - 還元 Fermi 準位 xi: 既定値 [-30, 10] を等間隔走査 - 判定: |(test - ref) / ref| <= tol 実装上の注意: - Seebeck_ZT.py 側は tau_pref を使うため、そこで得た equivalent l0 を Seebeck2.py 側へ渡して、同じ緩和時間スケールで比較する。 - 単位系の違いを吸収してから比較する。 * Seebeck_ZT.py: n[m^-3], S[V/K], sigma[S/m], mu[m^2/V/s] * Seebeck2.py : n[cm^-3], S[uV/K], sigma[S/cm], mu[cm^2/V/s] """ from __future__ import annotations import argparse import csv import importlib.util import math import sys from pathlib import Path from typing import Dict, List import numpy as np DEFAULT_R_LIST = [-0.5, 0.0, 0.5, 1.0, 1.5, 2.0] DEFAULT_XI_MIN = -30.0 DEFAULT_XI_MAX = 10.0 DEFAULT_NXI = 41 DEFAULT_T = 300.0 DEFAULT_M_EFF = 1.0 DEFAULT_MU_REF = 10.0 # cm^2/V/s DEFAULT_XI_REF = 0.0 DEFAULT_TOL = 1.0e-3 def parse_r_list(text: str) -> List[float]: vals = [] for token in text.split(","): token = token.strip() if token: vals.append(float(token)) if not vals: raise ValueError("r_list is empty") return vals def load_module_from_path(path: Path, module_name: str): print(f"\nLoad model: {str(path)}") spec = importlib.util.spec_from_file_location(module_name, str(path)) if spec is None or spec.loader is None: raise ImportError(f"Cannot load module from: {path}") module = importlib.util.module_from_spec(spec) spec.loader.exec_module(module) print(f" Fermi integral type: {module.fermi_integral_type}") return module def safe_rel_err(ref: float, val: float) -> float: if not math.isfinite(ref) or not math.isfinite(val): return math.inf if ref == 0.0: return 0.0 if val == 0.0 else math.inf return abs((val - ref) / ref) def configure_seebeck2(mod, T: float, m_eff: float, r: float, l0: float, carrier: str) -> None: """Seebeck2 側のグローバル状態を比較用に設定する。""" mod.T0 = T mod.klatt = 5.0 # 本比較では未使用だが、関数引数に必要 mod.dos.meeff = m_eff mod.dos.EC = 0.0 mod.dos.NC = mod.meff2NC_FEA(mod.dos.meeff, T) mod.dos.DC0 = mod.meff2DC0_FEA(mod.dos.meeff, T) mod.mobility.debug = 0 mod.mobility.use_simple = 0 mod.mobility.meff = m_eff mod.mobility.rfac = r mod.mobility.l0 = l0 mod.mobility.charge = -1.0 if carrier.lower().startswith("e") else 1.0 mod.mobility.set_scattering_parameters(mod.mobility.l0, mod.mobility.rfac) def calc_from_seebeck_zt(mod, xi: float, r: float, T: float, tau_pref: float, m_eff: float, carrier: str) -> Dict[str, float]: n_m3 = mod.electron_density_from_xi(xi, T=T, m_eff=m_eff) sigma_Sm = mod.sigma_from_xi_tau_pref(xi, r, T, tau_pref, m_eff=m_eff) mu_m2_Vs = mod.mobility_from_sigma_n(sigma_Sm, n_m3) S_V_K = mod.seebeck_from_xi_transport(xi, r, carrier=carrier) L_Wohm_K2 = mod.lorenz_from_xi_transport(xi, r) return { "ne_cm3": n_m3 / 1.0e6, "S_uVK": S_V_K * 1.0e6, "L": L_Wohm_K2, "mu_cm2_Vs": mu_m2_Vs * 1.0e4, "sigma_Scm": sigma_Sm / 100.0, } def calc_from_seebeck2(mod, xi: float, T: float) -> Dict[str, float]: EF_eV = xi * mod.kB * T / mod.e sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf = mod.dos.cal_transport_properteis( T, EF_eV, mod.mobility, mod.klatt, validate_error_str="compare_seebeck_results", ) return { "ne_cm3": n, "S_uVK": S, "L": L, "mu_cm2_Vs": mu, "sigma_Scm": sigma, } def summarize_metric(rows: List[dict], metric: str, tol: float) -> dict: worst = max(rows, key=lambda row: row[f"relerr_{metric}"]) max_relerr = worst[f"relerr_{metric}"] return { "metric": metric, "max_relerr": max_relerr, "worst_r": worst["r"], "worst_xi": worst["xi"], "ref": worst[f"zt_{metric}"], "test": worst[f"s2_{metric}"], "pass": max_relerr <= tol, } def write_csv(path: Path, fieldnames: List[str], rows: List[dict]) -> None: with path.open("w", newline="", encoding="utf-8-sig") as fp: writer = csv.DictWriter(fp, fieldnames=fieldnames) writer.writeheader() writer.writerows(rows) def main() -> int: parser = argparse.ArgumentParser( description="Compare Seebeck2.py and Seebeck_ZT.py over r and xi." ) parser.add_argument("--seebeck2", type=str, default="Seebeck2.py", help="Path to Seebeck2.py") parser.add_argument("--seebeck_zt", type=str, default="Seebeck_ZT.py", help="Path to Seebeck_ZT.py") parser.add_argument("--T", type=float, default=DEFAULT_T, help="Temperature [K]") parser.add_argument("--m_eff", type=float, default=DEFAULT_M_EFF, help="Effective mass m*/m0") parser.add_argument("--mu_ref", type=float, default=DEFAULT_MU_REF, help="Reference mobility [cm^2/V/s]") parser.add_argument("--xi_ref", type=float, default=DEFAULT_XI_REF, help="Reference reduced Fermi level") parser.add_argument("--carrier", type=str, default="electron", choices=["electron", "hole"], help="Carrier type") parser.add_argument( "--r_list", type=str, default=",".join(str(x) for x in DEFAULT_R_LIST), help="Comma-separated r list", ) parser.add_argument("--xi_min", type=float, default=DEFAULT_XI_MIN, help="Minimum xi") parser.add_argument("--xi_max", type=float, default=DEFAULT_XI_MAX, help="Maximum xi") parser.add_argument("--nxi", type=int, default=DEFAULT_NXI, help="Number of xi points") parser.add_argument("--tol", type=float, default=DEFAULT_TOL, help="Tolerance for absolute relative error") parser.add_argument("--detail_csv", type=str, default="compare_seebeck_detail.csv", help="Detailed CSV output") parser.add_argument("--summary_csv", type=str, default="compare_seebeck_summary.csv", help="Summary CSV output") args = parser.parse_args() s2_path = Path(args.seebeck2) zt_path = Path(args.seebeck_zt) if not s2_path.exists(): print(f"ERROR: Seebeck2 file not found: {s2_path}", file=sys.stderr) return 2 if not zt_path.exists(): print(f"ERROR: Seebeck_ZT file not found: {zt_path}", file=sys.stderr) return 2 try: seebeck2 = load_module_from_path(s2_path, "seebeck2_compare_mod") except ModuleNotFoundError as exc: print( "ERROR: Failed to import Seebeck2.py.\n" f"Missing module: {exc.name}\n" "Seebeck2.py depends on your tklib environment, so run this script in the same environment where Seebeck2.py works.", file=sys.stderr, ) return 2 try: seebeck_zt = load_module_from_path(zt_path, "seebeck_zt_compare_mod") except ModuleNotFoundError as exc: print( "ERROR: Failed to import Seebeck_ZT.py.\n" f"Missing module: {exc.name}", file=sys.stderr, ) return 2 r_list = parse_r_list(args.r_list) xi_values = np.linspace(args.xi_min, args.xi_max, args.nxi) print("Comparison settings") print(f" Seebeck2 : {s2_path}") print(f" Seebeck_ZT : {zt_path}") print(f" T : {args.T}") print(f" m_eff : {args.m_eff}") print(f" mu_ref : {args.mu_ref} cm^2/V/s") print(f" xi_ref : {args.xi_ref}") print(f" carrier : {args.carrier}") print(f" r_list : {r_list}") print(f" xi range : [{args.xi_min}, {args.xi_max}] n={args.nxi}") print(f" tol : {args.tol}") print() detail_rows: List[dict] = [] summary_rows: List[dict] = [] metrics = ["ne_cm3", "S_uVK", "L", "mu_cm2_Vs"] overall_pass = True for r in r_list: tau_pref = seebeck_zt.tau_pref_from_mu_ref(args.mu_ref, args.xi_ref, r, args.T, m_eff=args.m_eff) l0_equiv = seebeck_zt.equivalent_l0_from_tau_pref(tau_pref, r, m_eff=args.m_eff) configure_seebeck2(seebeck2, args.T, args.m_eff, r, l0_equiv, args.carrier) rows_this_r = [] for xi in xi_values: ref = calc_from_seebeck_zt(seebeck_zt, float(xi), r, args.T, tau_pref, args.m_eff, args.carrier) test = calc_from_seebeck2(seebeck2, float(xi), args.T) row = { "r": r, "xi": float(xi), "tau_pref": tau_pref, "l0_equiv_m": l0_equiv, } for key in ["ne_cm3", "S_uVK", "L", "mu_cm2_Vs", "sigma_Scm"]: row[f"zt_{key}"] = ref[key] row[f"s2_{key}"] = test[key] for key in metrics: row[f"relerr_{key}"] = safe_rel_err(ref[key], test[key]) row["pass"] = all(row[f"relerr_{key}"] <= args.tol for key in metrics) if not row["pass"]: overall_pass = False detail_rows.append(row) rows_this_r.append(row) print(f"r = {r:g}") print(f" tau_pref = {tau_pref:.8e}") print(f" l0_equiv = {l0_equiv:.8e} m") print(" Detailed comparison") for row in rows_this_r: print(f" xi = {row['xi']: .6f}") for key in metrics: ok_ng = "PASS" if row[f"relerr_{key}"] <= args.tol else "FAIL" print( f" {key:10s} " f"Seebeck2={row[f's2_{key}']:+.12e} " f"Seebeck_ZT={row[f'zt_{key}']:+.12e} " f"relerr={row[f'relerr_{key}']:+.6e} [{ok_ng}]" ) print(" Summary for max relative error in this r") for key in metrics: worst = max(rows_this_r, key=lambda row: row[f"relerr_{key}"]) ok_ng = "PASS" if worst[f"relerr_{key}"] <= args.tol else "FAIL" print( f" {key:10s}: max_relerr = {worst[f'relerr_{key}']:.6e} " f"at xi = {worst['xi']: .6f} " f"Seebeck2={worst[f's2_{key}']:+.12e} " f"Seebeck_ZT={worst[f'zt_{key}']:+.12e} [{ok_ng}]" ) summary_rows.append( { "r": r, "metric": key, "tau_pref": tau_pref, "l0_equiv_m": l0_equiv, "max_relerr": worst[f"relerr_{key}"], "worst_xi": worst["xi"], "ref_value": worst[f"zt_{key}"], "test_value": worst[f"s2_{key}"], "pass": worst[f"relerr_{key}"] <= args.tol, "tol": args.tol, } ) print() detail_csv = Path(args.detail_csv) summary_csv = Path(args.summary_csv) write_csv( detail_csv, [ "r", "xi", "tau_pref", "l0_equiv_m", "zt_ne_cm3", "s2_ne_cm3", "relerr_ne_cm3", "zt_S_uVK", "s2_S_uVK", "relerr_S_uVK", "zt_L", "s2_L", "relerr_L", "zt_mu_cm2_Vs", "s2_mu_cm2_Vs", "relerr_mu_cm2_Vs", "zt_sigma_Scm", "s2_sigma_Scm", "pass", ], detail_rows, ) write_csv( summary_csv, [ "r", "metric", "tau_pref", "l0_equiv_m", "max_relerr", "worst_xi", "ref_value", "test_value", "pass", "tol", ], summary_rows, ) metric_summaries = [summarize_metric(detail_rows, metric, args.tol) for metric in metrics] overall_worst = max( metric_summaries, key=lambda item: item["max_relerr"], ) print("Final summary by maximum relative error") for item in metric_summaries: ok_ng = "PASS" if item["pass"] else "FAIL" print( f" {item['metric']:10s}: max_relerr = {item['max_relerr']:.6e} " f"at r = {item['worst_r']:g}, xi = {item['worst_xi']: .6f} " f"Seebeck2={item['test']:+.12e} Seebeck_ZT={item['ref']:+.12e} [{ok_ng}]" ) print() print("Overall worst case") print(f" metric = {overall_worst['metric']}") print(f" max_relerr = {overall_worst['max_relerr']:.6e}") print(f" r = {overall_worst['worst_r']:g}") print(f" xi = {overall_worst['worst_xi']: .6f}") print(f" Seebeck2 = {overall_worst['test']:.12e}") print(f" Seebeck_ZT = {overall_worst['ref']:.12e}") print(f" result = {'PASS' if overall_worst['pass'] else 'FAIL'}") print() print(f"Detailed results saved to: {detail_csv}") print(f"Summary results saved to: {summary_csv}") print() print("OVERALL:", "PASS" if overall_pass else "FAIL") return 0 if overall_pass else 1 if __name__ == "__main__": raise SystemExit(main())