Seebeck_ZT.py ダウンロード/コピー
Seebeck_ZT.py
Seebeck_ZT.py
1"""
2熱電輸送特性計算モジュール
3
4概要:
5ゼーベック係数、ローレンツ因子、電気伝導度、移動度、電子熱伝導度、パワーファクター、ZT因子を計算し、
6キャリア濃度依存性をプロット・Excelファイルに保存します。
7
8詳細説明:
9このモジュールは、3次元放物線バンドモデルと緩和時間近似（tau(E) ∝ E^(r-1/2)）に基づき、
10フェルミ・ディラック積分を用いた厳密な輸送積分によって熱電輸送特性を評価します。
11コマンドライン引数を通じて、温度、有効質量、参照移動度、格子熱伝導率、キャリア種別、散乱因子rのリストなどを設定できます。
12
13関連リンク:
14Seebeck_ZT_usage
15"""
16import numpy as np
17import matplotlib.pyplot as plt
18from scipy.integrate import quad
19from scipy.special import expit, gamma
20from functools import lru_cache
21import argparse
22from openpyxl import Workbook
23
24# =========================
25# constants (SI)
26# =========================
27kB = 1.380649e-23 # J/K
28e = 1.602176634e-19 # C
29h = 6.62607015e-34 # J s
30hbar = h / (2.0 * np.pi)
31m0 = 9.1093837015e-31 # kg
32
33r"""
34散乱因子 r は tau(E) ∝ E^(r-1/2) で定義
35"""
36
37figsize = (10, 10)
38
39# ============================================================
40# Fermi–Dirac integral (NO 1/Gamma(j+1))
41# Fj(j, xi) = ∫_0^∞ x^j / (1+exp(x-xi)) dx
42# ============================================================
43@lru_cache(maxsize=40000)
44def Fj(j, xi):
45 """
46 概要:
47 フェルミ・ディラック積分 Fj(j, xi) = ∫_0^∞ x^j / (1+exp(x-xi)) dx を計算します。
48
49 詳細説明:
50 SciPyのquad関数を使用して数値積分を実行します。結果はlru_cacheでキャッシュされ、
51 同じ引数での再計算を回避し、パフォーマンスを向上させます。
52 この関数は、ガンマ関数による正規化（1/Gamma(j+1)）は含まれていません。
53
54 引数:
55 :param j: 積分の次数。
56 :type j: float
57 :param xi: 還元フェルミ準位 (reduced Fermi level)。
58 :type xi: float
59 戻り値:
60 :returns: フェルミ・ディラック積分 Fj(j, xi) の値。
61 :rtype: float
62 """
63 j = float(j)
64 xi = float(np.round(xi, 12))
65 integrand = lambda x: (x**j) * expit(xi - x)
66 val, _ = quad(integrand, 0.0, np.inf, epsabs=1e-10, epsrel=1e-10, limit=400)
67 return val
68
69
70def Nc_3D(T, m_eff=1.0):
71 """
72 概要:
73 3次元放物線バンドにおける有効状態密度 Nc を計算します。
74
75 詳細説明:
76 以下の公式を用いて有効状態密度を計算します。
77 Nc = 2 * (2π m* kB T / h^2)^(3/2)
78 ここで、m* は有効質量、kB はボルツマン定数、T は温度、h はプランク定数です。
79
80 引数:
81 :param T: 温度 [K]。
82 :type T: float
83 :param m_eff: 有効質量の比率 (m*/m0)。m0 は自由電子質量です。デフォルトは1.0。
84 :type m_eff: float, optional
85 戻り値:
86 :returns: 有効状態密度 Nc [m^-3]。
87 :rtype: float
88 """
89 mstar = m_eff * m0
90 return 2.0 * (2.0 * np.pi * mstar * kB * T / (h**2))**1.5
91
92
93def electron_density_from_xi(xi, T, m_eff=1.0):
94 """
95 概要:
96 還元フェルミ準位 xi から電子濃度 n を計算します。
97
98 詳細説明:
99 電子濃度は以下の式で計算されます。
100 n [m^-3] = Nc * F^{std}_{1/2}(xi)
101 ここで、Nc は有効状態密度、F^{std}_{1/2}(xi) = Fj(1/2, xi) / Gamma(3/2) は
102 標準化されたフェルミ・ディラック積分です。
103
104 引数:
105 :param xi: 還元フェルミ準位。
106 :type xi: float
107 :param T: 温度 [K]。
108 :type T: float
109 :param m_eff: 有効質量の比率 (m*/m0)。デフォルトは1.0。
110 :type m_eff: float, optional
111 戻り値:
112 :returns: 電子濃度 n [m^-3]。
113 :rtype: float
114 """
115 Nc = Nc_3D(T, m_eff=m_eff)
116 F12_std = Fj(0.5, xi) / gamma(1.5)
117 return Nc * F12_std
118
119
120# ============================================================
121# Transport integrals for tau(E) ∝ E^(r-1/2)
122# I0 = (r+1) F_r
123# I1 = (r+2) F_{r+1} - xi (r+1) F_r
124# I2 = (r+3) F_{r+2} - 2 xi (r+2) F_{r+1} + xi^2 (r+1) F_r
125# ============================================================
126def transport_I0(xi, r):
127 """
128 概要:
129 輸送積分 I0 を計算します。
130
131 詳細説明:
132 tau(E) ∝ E^(r-1/2) の緩和時間に対する輸送積分の一つである I0 を計算します。
133 I0 = (r + 1) * F_r(xi) の式で定義されます。
134
135 引数:
136 :param xi: 還元フェルミ準位。
137 :type xi: float
138 :param r: 散乱因子。緩和時間のエネルギー依存性 E^(r-1/2) を定義します。
139 :type r: float
140 戻り値:
141 :returns: 輸送積分 I0 の値。
142 :rtype: float
143 """
144 return (r + 1.0) * Fj(r, xi)
145
146
147def transport_I1(xi, r):
148 """
149 概要:
150 輸送積分 I1 を計算します。
151
152 詳細説明:
153 tau(E) ∝ E^(r-1/2) の緩和時間に対する輸送積分の一つである I1 を計算します。
154 I1 = (r + 2) * F_{r+1}(xi) - xi * (r + 1) * F_r(xi) の式で定義されます。
155
156 引数:
157 :param xi: 還元フェルミ準位。
158 :type xi: float
159 :param r: 散乱因子。緩和時間のエネルギー依存性 E^(r-1/2) を定義します。
160 :type r: float
161 戻り値:
162 :returns: 輸送積分 I1 の値。
163 :rtype: float
164 """
165 return (r + 2.0) * Fj(r + 1.0, xi) - xi * (r + 1.0) * Fj(r, xi)
166
167
168def transport_I2(xi, r):
169 """
170 概要:
171 輸送積分 I2 を計算します。
172
173 詳細説明:
174 tau(E) ∝ E^(r-1/2) の緩和時間に対する輸送積分の一つである I2 を計算します。
175 I2 = (r + 3) * F_{r+2}(xi) - 2 * xi * (r + 2) * F_{r+1}(xi) + xi^2 * (r + 1) * F_r(xi)
176 の式で定義されます。
177
178 引数:
179 :param xi: 還元フェルミ準位。
180 :type xi: float
181 :param r: 散乱因子。緩和時間のエネルギー依存性 E^(r-1/2) を定義します。
182 :type r: float
183 戻り値:
184 :returns: 輸送積分 I2 の値。
185 :rtype: float
186 """
187 return (
188 (r + 3.0) * Fj(r + 2.0, xi)
189 - 2.0 * xi * (r + 2.0) * Fj(r + 1.0, xi)
190 + (xi**2) * (r + 1.0) * Fj(r, xi)
191 )
192
193
194def transport_prefactor_sigma(T, m_eff=1.0):
195 """
196 概要:
197 電気伝導度 sigma の計算における前因子 A_sigma を計算します。
198
199 詳細説明:
200 電気伝導度は sigma = A_sigma * tau_pref * (kB T)^(r+1) * I0 の形で表されます。
201 この関数は、tau(E) = tau_pref * E^(r-1/2) と定義される緩和時間において、
202 以下の物理定数を含む前因子 A_sigma を計算します。
203 A_sigma [SI] = e^2 * 2^(3/2) * sqrt(m*) / (3 π^2 ħ^3)
204 ここで、e は素電荷、m* は有効質量、ħ は換算プランク定数です。
205
206 引数:
207 :param T: 温度 [K]。
208 :type T: float
209 :param m_eff: 有効質量の比率 (m*/m0)。m0 は自由電子質量です。デフォルトは1.0。
210 :type m_eff: float, optional
211 戻り値:
212 :returns: 電気伝導度の前因子 A_sigma [SI単位]。
213 :rtype: float
214 """
215 mstar = m_eff * m0
216 return e**2 * (2.0**1.5) * np.sqrt(mstar) / (3.0 * np.pi**2 * hbar**3)
217
218
219def sigma_from_xi_tau_pref(xi, r, T, tau_pref, m_eff=1.0):
220 """
221 概要:
222 還元フェルミ準位 xi、散乱因子 r、温度 T、緩和時間の前因子 tau_pref から
223 電気伝導度 sigma を計算します。
224
225 詳細説明:
226 電気伝導度は以下の式で計算されます。
227 sigma = A_sigma * tau_pref * (kB T)^(r + 1.0) * I0
228 ここで、A_sigma は transport_prefactor_sigma で計算される前因子、
229 I0 は transport_I0 で計算される輸送積分です。
230
231 引数:
232 :param xi: 還元フェルミ準位。
233 :type xi: float
234 :param r: 散乱因子。
235 :type r: float
236 :param T: 温度 [K]。
237 :type T: float
238 :param tau_pref: 緩和時間の前因子 tau_pref [s / J^(r-1/2)]。
239 :type tau_pref: float
240 :param m_eff: 有効質量の比率 (m*/m0)。m0 は自由電子質量です。デフォルトは1.0。
241 :type m_eff: float, optional
242 戻り値:
243 :returns: 電気伝導度 sigma [S/m]。
244 :rtype: float
245 """
246 A = transport_prefactor_sigma(T, m_eff=m_eff)
247 I0 = transport_I0(xi, r)
248 return A * tau_pref * (kB * T)**(r + 1.0) * I0
249
250
251def seebeck_from_xi_transport(xi, r, carrier="electron"):
252 """
253 概要:
254 還元フェルミ準位 xi と散乱因子 r からゼーベック係数 S を計算します。
255
256 詳細説明:
257 ゼーベック係数は厳密な輸送積分を用いて以下の式で計算されます。
258 S = (kB/q) * I1/I0
259 q はキャリアの種類によって異なり、電子の場合は -e、正孔の場合は +e となります。
260
261 引数:
262 :param xi: 還元フェルミ準位。
263 :type xi: float
264 :param r: 散乱因子。tau(E) ∝ E^(r-1/2) で定義されます。
265 :type r: float
266 :param carrier: キャリアの種類。"electron"または"hole"を指定します。デフォルトは"electron"。
267 :type carrier: str, optional
268 戻り値:
269 :returns: ゼーベック係数 S [V/K]。
270 :rtype: float
271 例外:
272 :raises ValueError: carrier が "electron" または "hole" 以外の場合。
273 """
274 if carrier.lower().startswith("e"):
275 q = -e
276 elif carrier.lower().startswith("h"):
277 q = +e
278 else:
279 raise ValueError("carrier must be 'electron' or 'hole'")
280
281 I0 = transport_I0(xi, r)
282 I1 = transport_I1(xi, r)
283 return (kB / q) * (I1 / I0)
284
285
286def lorenz_from_xi_transport(xi, r):
287 """
288 概要:
289 還元フェルミ準位 xi と散乱因子 r からローレンツ因子 L を計算します。
290
291 詳細説明:
292 ローレンツ因子は厳密な輸送積分を用いて以下の式で計算されます。
293 L = (kB/e)^2 * ( I2/I0 - (I1/I0)^2 )
294
295 引数:
296 :param xi: 還元フェルミ準位。
297 :type xi: float
298 :param r: 散乱因子。
299 :type r: float
300 戻り値:
301 :returns: ローレンツ因子 L [V^2/K^2 (または W ohm / K^2)]。
302 :rtype: float
303 """
304 I0 = transport_I0(xi, r)
305 I1 = transport_I1(xi, r)
306 I2 = transport_I2(xi, r)
307 return (kB / e)**2 * (I2 / I0 - (I1 / I0)**2)
308
309
310def mobility_from_sigma_n(sigma, n):
311 """
312 概要:
313 電気伝導度 sigma と電子濃度 n から移動度 mu を計算します。
314
315 詳細説明:
316 電気伝導度 sigma、電子濃度 n、素電荷 e、移動度 mu の間に成り立つ関係式
317 sigma = n * e * mu を用いて移動度を計算します。
318 したがって、mu = sigma / (n * e) となります。
319
320 引数:
321 :param sigma: 電気伝導度 [S/m]。
322 :type sigma: float
323 :param n: 電子濃度 [m^-3]。
324 :type n: float
325 戻り値:
326 :returns: 移動度 mu [m^2/V/s]。
327 :rtype: float
328 """
329 return sigma / (n * e)
330
331
332def tau_pref_from_mu_ref(mu_ref_cm2_Vs, xi_ref, r, T, m_eff=1.0):
333 """
334 概要:
335 参照移動度 mu_ref と参照還元フェルミ準位 xi_ref に基づいて、
336 緩和時間の前因子 tau_pref を決定します。
337
338 詳細説明:
339 緩和時間 tau(E) が tau(E) = tau_pref * E^(r-1/2) で定義されるとき、
340 指定された xi_ref と T における移動度が mu_ref となるように
341 tau_pref を逆算します。
342 参照移動度 mu_ref_cm2_Vs は cm^2/V/s 単位で与えられますが、計算ではSI単位に変換されます。
343
344 引数:
345 :param mu_ref_cm2_Vs: 参照移動度 [cm^2/V/s]。
346 :type mu_ref_cm2_Vs: float
347 :param xi_ref: 参照還元フェルミ準位。この xi で参照移動度が再現されます。
348 :type xi_ref: float
349 :param r: 散乱因子。
350 :type r: float
351 :param T: 温度 [K]。
352 :type T: float
353 :param m_eff: 有効質量の比率 (m*/m0)。m0 は自由電子質量です。デフォルトは1.0。
354 :type m_eff: float, optional
355 戻り値:
356 :returns: 緩和時間の前因子 tau_pref [s / J^(r-1/2)]。
357 :rtype: float
358 """
359 mu_ref_SI = mu_ref_cm2_Vs * 1e-4 # cm^2/V/s -> m^2/V/s
360 n_ref = electron_density_from_xi(xi_ref, T=T, m_eff=m_eff)
361 sigma_ref = n_ref * e * mu_ref_SI
362
363 A = transport_prefactor_sigma(T, m_eff=m_eff)
364 I0_ref = transport_I0(xi_ref, r)
365
366 tau_pref = sigma_ref / (A * (kB * T)**(r + 1.0) * I0_ref)
367 return tau_pref
368
369
370def tau_at_energy(E_J, tau_pref, r):
371 """
372 概要:
373 特定のエネルギー E_J における緩和時間 tau(E_J) を計算します。
374
375 詳細説明:
376 緩和時間は以下の式で定義されます。
377 tau(E) = tau_pref * E^(r-1/2)
378 E_J が0以下の場合は np.nan を返します。
379
380 引数:
381 :param E_J: エネルギー [J]。
382 :type E_J: float
383 :param tau_pref: 緩和時間の前因子 [s / J^(r-1/2)]。
384 :type tau_pref: float
385 :param r: 散乱因子。
386 :type r: float
387 戻り値:
388 :returns: エネルギー E_J における緩和時間 tau(E_J) [s]。
389 :rtype: float
390 """
391 if E_J <= 0.0:
392 return np.nan
393 return tau_pref * (E_J**(r - 0.5))
394
395
396def equivalent_l0_from_tau_pref(tau_pref, r, m_eff=1.0):
397 """
398 概要:
399 緩和時間の前因子 tau_pref から等価な平均自由行程 l0 を計算します。
400
401 詳細説明:
402 一部のスライドでの慣例的な定義
403 tau(E) = sqrt(m*/2) * [ l0(T) * e_chg^(-r) ] * E^(r-1/2)
404 から、l0(T) を逆算します。
405 したがって、l0(T) = tau_pref * sqrt(2/m*) * e_chg^r となります。
406
407 引数:
408 :param tau_pref: 緩和時間の前因子 [s / J^(r-1/2)]。
409 :type tau_pref: float
410 :param r: 散乱因子。
411 :type r: float
412 :param m_eff: 有効質量の比率 (m*/m0)。m0 は自由電子質量です。デフォルトは1.0。
413 :type m_eff: float, optional
414 戻り値:
415 :returns: 等価な平均自由行程 l0 [m]。
416 :rtype: float
417 """
418 mstar = m_eff * m0
419 return tau_pref * np.sqrt(2.0 / mstar) * (e**r)
420
421
422def save_results_to_excel(outfile, meta_rows, data_rows):
423 """
424 概要:
425 計算結果とメタデータをExcelファイルに保存します。
426
427 詳細説明:
428 指定された出力ファイル名で新しいExcelワークブックを作成し、
429 "metadata"と"transport_vs_Ne"の2つのシートにデータを書き込みます。
430 "metadata"シートには計算条件などのメタデータが保存され、
431 "transport_vs_Ne"シートにはキャリア濃度に対する輸送特性の計算結果が保存されます。
432
433 引数:
434 :param outfile: 出力するExcelファイル名。例: "results.xlsx"。
435 :type outfile: str
436 :param meta_rows: メタデータを含む行のリスト。各要素はExcelの1行に対応するリストです。
437 :type meta_rows: list
438 :param data_rows: 計算結果データを含む行のリスト。各要素はExcelの1行に対応するリストです。
439 :type data_rows: list
440 戻り値:
441 :returns: なし
442 :rtype: None
443 """
444 wb = Workbook()
445 ws_meta = wb.active
446 ws_meta.title = "metadata"
447
448 for row in meta_rows:
449 ws_meta.append(row)
450
451 ws = wb.create_sheet("transport_vs_Ne")
452 ws.append([
453 "r",
454 "xi",
455 "tau_pref_s_per_J^(r-1/2)",
456 "tau_at_kBT_s",
457 "l0_equiv_m",
458 "Ne_m^-3",
459 "Ne_cm^-3",
460 "S_V_per_K",
461 "S_uV_per_K",
462 "L_Wohm_per_K2",
463 "sigma_S_per_m",
464 "mu_m2_per_Vs",
465 "mu_cm2_per_Vs",
466 "kappa_e_W_per_mK",
467 "PF_W_per_mK2",
468 "PF_mW_per_mK2",
469 "ZT",
470 ])
471 for row in data_rows:
472 ws.append(row)
473
474 wb.save(outfile)
475
476
477def main():
478 """
479 概要:
480 スクリプトの主要な実行フローを定義します。コマンドライン引数を解析し、
481 熱電輸送特性を計算・プロット・保存します。
482
483 詳細説明:
484 argparse モジュールを使用してコマンドライン引数（温度、有効質量、参照移動度、
485 格子熱伝導率、キャリア種別、還元フェルミ準位の範囲と点数、散乱因子rのリスト、
486 出力ファイル名など）を解析します。
487 次に、指定されたパラメータに基づいて、各散乱因子 r に対して輸送特性（電子濃度、
488 ゼーベック係数、ローレンツ因子、電気伝導度、移動度、電子熱伝導度、パワーファクター、ZT因子）を計算します。
489 計算されたデータは save_results_to_excel 関数を使用してExcelファイルに保存されます。
490 最後に、計算結果は matplotlib を用いてキャリア濃度に対する各特性のグラフとしてプロットされ、表示されます。
491
492 戻り値:
493 :returns: なし
494 :rtype: None
495 """
496 parser = argparse.ArgumentParser(
497 description=(
498 "S, L, sigma, mu, kappa_e, PF, ZT vs Ne using exact transport integrals "
499 "for tau(E) ∝ E^(r-1/2) in a 3D parabolic band."
500 )
501 )
502 parser.add_argument("--T", type=float, default=300.0, help="Temperature [K]")
503 parser.add_argument("--m_eff", type=float, default=1.0, help="Effective mass m*/m0")
504 parser.add_argument("--mu_ref", type=float, default=10.0,
505 help="Reference mobility at xi_ref [cm^2/V/s]")
506 parser.add_argument("--xi_ref", type=float, default=0.0,
507 help="Reference reduced Fermi level xi where mu_ref is reproduced")
508 parser.add_argument("--klatt", type=float, default=5.0,
509 help="Lattice thermal conductivity [W/m/K] (default: 5)")
510 parser.add_argument("--carrier", type=str, default="electron",
511 choices=["electron", "hole"], help="Carrier type")
512 parser.add_argument("--xi_min", type=float, default=-5.0, help="Minimum xi")
513 parser.add_argument("--xi_max", type=float, default=40.0, help="Maximum xi")
514 parser.add_argument("--nxi", type=int, default=181, help="Number of xi points")
515 parser.add_argument("--r_list", type=str, default="2.0,1.5,1.0,0.5,0.0,-0.5",
516 help="Comma-separated r list")
517 parser.add_argument("--outfile", type=str, default="transport_vs_Ne.xlsx",
518 help="Excel output filename")
519 args = parser.parse_args()
520
521 T = args.T
522 m_eff = args.m_eff
523 mu_ref_cm2_Vs = args.mu_ref
524 xi_ref = args.xi_ref
525 klatt = args.klatt
526 carrier = args.carrier
527 r_list = [float(x.strip()) for x in args.r_list.split(",")]
528 xi_range = np.linspace(args.xi_min, args.xi_max, args.nxi)
529 outfile = args.outfile
530
531 print("Convention: tau(E) ∝ E^(r-1/2)")
532 print(f"T = {T:g} K")
533 print(f"m_eff = {m_eff:g} m0")
534 print(f"mu_ref = {mu_ref_cm2_Vs:g} cm^2/V/s")
535 print(f"xi_ref = {xi_ref:g}")
536 print(f"klatt = {klatt:g} W/m/K")
537 print(f"carrier = {carrier}")
538 print(f"outfile = {outfile}")
539 print()
540
541 all_results = []
542
543 for r in r_list:
544 tau_pref = tau_pref_from_mu_ref(mu_ref_cm2_Vs, xi_ref, r, T, m_eff=m_eff)
545 tau_kBT = tau_at_energy(kB * T, tau_pref, r)
546 l0_equiv = equivalent_l0_from_tau_pref(tau_pref, r, m_eff=m_eff)
547
548 print(f"r = {r:g}")
549 print(f" tau_pref = {tau_pref:.6e} [s / J^({r:g}-1/2)]")
550 print(f" tau(E=kBT) = {tau_kBT:.6e} [s]")
551 print(f" equivalent l0 = {l0_equiv:.6e} [m] (slide convention)")
552
553 n_list = []
554 S_list = []
555 L_list = []
556 sig_list = []
557 mu_list = []
558 ke_list = []
559 PF_list = []
560 ZT_list = []
561
562 for xi in xi_range:
563 n = electron_density_from_xi(xi, T=T, m_eff=m_eff) # [m^-3]
564 sigma = sigma_from_xi_tau_pref(xi, r, T, tau_pref, m_eff=m_eff) # [S/m]
565 mu_SI = mobility_from_sigma_n(sigma, n) # [m^2/V/s]
566 mu_cm2_Vs = mu_SI * 1e4
567
568 S = seebeck_from_xi_transport(xi, r, carrier=carrier) # [V/K]
569 L = lorenz_from_xi_transport(xi, r) # [V^2/K^2]
570 kappa_e = L * sigma * T # [W/m/K]
571 PF = S**2 * sigma # [W/m/K^2]
572 ZT = PF * T / (kappa_e + klatt)
573
574 n_list.append(n)
575 S_list.append(S)
576 L_list.append(L)
577 sig_list.append(sigma)
578 mu_list.append(mu_cm2_Vs)
579 ke_list.append(kappa_e)
580 PF_list.append(PF)
581 ZT_list.append(ZT)
582
583 all_results.append([
584 r,
585 xi,
586 tau_pref,
587 tau_kBT,
588 l0_equiv,
589 n,
590 n / 1e6,
591 S,
592 S * 1e6,
593 L,
594 sigma,
595 mu_SI,
596 mu_cm2_Vs,
597 kappa_e,
598 PF,
599 PF * 1e3,
600 ZT,
601 ])
602
603 # 配列化して保持
604 # この行は、ループ内でリストに追加された要素を再度リストに代入するもので、
605 # 結果的にはall_resultsのリスト構造に影響を与えません。
606 # 論理的な意味は薄いですが、既存コードの変更禁止ルールに従い保持します。
607 all_results[-len(xi_range):] = all_results[-len(xi_range):]
608
609 # ---------- save Excel before plotting ----------
610 meta_rows = [
611 ["Convention", "tau(E) ∝ E^(r-1/2)"],
612 ["T_K", T],
613 ["m_eff_over_m0", m_eff],
614 ["mu_ref_cm2_per_Vs", mu_ref_cm2_Vs],
615 ["xi_ref", xi_ref],
616 ["klatt_W_per_mK", klatt],
617 ["carrier", carrier],
618 ["r_list", ",".join(str(r) for r in r_list)],
619 ["xi_min", args.xi_min],
620 ["xi_max", args.xi_max],
621 ["nxi", args.nxi],
622 ]
623 save_results_to_excel(outfile, meta_rows, all_results)
624 print(f"\nSaved calculation results to [{outfile}]\n")
625
626 # ---------- plot ----------
627 fig, axes = plt.subplots(4, 2, figsize=figsize)
628 axS = axes[0, 0]
629 axL = axes[0, 1]
630 axSig = axes[1, 0]
631 axMu = axes[1, 1]
632 axKe = axes[2, 0]
633 axPF = axes[2, 1]
634 axZT = axes[3, 0]
635 axBlank = axes[3, 1]
636
637 for r in r_list:
638 tau_pref = tau_pref_from_mu_ref(mu_ref_cm2_Vs, xi_ref, r, T, m_eff=m_eff)
639
640 n_list = []
641 S_list = []
642 L_list = []
643 sig_list = []
644 mu_list = []
645 ke_list = []
646 PF_list = []
647 ZT_list = []
648
649 for xi in xi_range:
650 n = electron_density_from_xi(xi, T=T, m_eff=m_eff)
651 sigma = sigma_from_xi_tau_pref(xi, r, T, tau_pref, m_eff=m_eff)
652 mu_SI = mobility_from_sigma_n(sigma, n)
653 S = seebeck_from_xi_transport(xi, r, carrier=carrier)
654 L = lorenz_from_xi_transport(xi, r)
655 kappa_e = L * sigma * T
656 PF = S**2 * sigma
657 ZT = PF * T / (kappa_e + klatt)
658
659 n_list.append(n)
660 S_list.append(S)
661 L_list.append(L)
662 sig_list.append(sigma)
663 mu_list.append(mu_SI * 1e4)
664 ke_list.append(kappa_e)
665 PF_list.append(PF)
666 ZT_list.append(ZT)
667
668 n_arr_cm3 = np.array(n_list) / 1e6
669 S_arr_uVK = np.array(S_list) * 1e6
670 L_arr = np.array(L_list)
671 sig_arr = np.array(sig_list)
672 mu_arr = np.array(mu_list)
673 ke_arr = np.array(ke_list)
674 PF_arr = np.array(PF_list) * 1e3
675 ZT_arr = np.array(ZT_list)
676
677 label = f"r = {r:g}"
678
679 axS.plot(n_arr_cm3, S_arr_uVK, label=label)
680 axL.plot(n_arr_cm3, L_arr, label=label)
681 axSig.plot(n_arr_cm3, sig_arr, label=label)
682 axMu.plot(n_arr_cm3, mu_arr, label=label)
683 axKe.plot(n_arr_cm3, ke_arr, label=label)
684 axPF.plot(n_arr_cm3, PF_arr, label=label)
685 axZT.plot(n_arr_cm3, ZT_arr, label=label)
686
687 for ax in [axS, axL, axSig, axMu, axKe, axPF, axZT]:
688 ax.set_xscale("log")
689 ax.grid(True, which="both", alpha=0.25)
690 ax.set_xlabel(r"$N_e$ (cm$^{-3}$)")
691
692 axS.axhline(0.0, linewidth=1.0, alpha=0.4)
693 axS.set_ylabel(r"$S$ ($\mu$V/K)")
694 axS.set_title("Seebeck coefficient")
695
696 axL.set_ylabel(r"$L$ (V$^2$/K$^2$)")
697 axL.set_title("Lorenz factor")
698
699 axSig.set_ylabel(r"$\sigma$ (S/m)")
700 axSig.set_title("Electrical conductivity")
701
702 axMu.set_ylabel(r"$\mu$ (cm$^2$/V/s)")
703 axMu.set_title(r"Mobility ($\tau(E)\propto E^{r-1/2}$)")
704
705 axKe.set_ylabel(r"$\kappa_e$ (W/m/K)")
706 axKe.set_title("Electronic thermal conductivity")
707
708 axPF.set_ylabel(r"$PF=S^2\sigma$ (mW/m/K$^2$)")
709 axPF.set_title("Power factor")
710
711 axZT.set_ylabel(r"$ZT$")
712 axZT.set_title(rf"$ZT$ ($\kappa_{{latt}}$={klatt:g} W/m/K)")
713
714 axBlank.axis("off")
715
716 handles, labels = axS.get_legend_handles_labels()
717 fig.legend(handles, labels, loc="upper center", ncol=min(len(r_list), 6))
718 fig.suptitle(
719 f"Transport vs Ne (T={T:g} K, m*/m0={m_eff:g}, "
720 f"mu_ref={mu_ref_cm2_Vs:g} cm$^2$/V/s at xi_ref={xi_ref:g}, "
721 r"$\tau(E)\propto E^{r-1/2}$)",
722 y=0.995
723 )
724 plt.tight_layout(rect=[0, 0, 1, 0.96])
725 plt.show()
726
727
728if __name__ == "__main__":
729 main()