"""
SPBモデルを用いたゼーベック係数からのローレンツ数と電子熱伝導率の推定スクリプト
このスクリプトは、Single Parabolic Band (SPB) モデルを用いて、与えられた電気伝導率と
ゼーベック係数からフェルミ準位を推定し、ローレンツ数および電子熱伝導率を計算します。
主に熱電材料の特性評価に利用されます。
関連リンク: :doc:`S2L_usage`
"""
#間違っていても責任はとりません、野元
#Single parabric band modelを使ってLorentz数、Keを計算します。
# -*- coding: utf-8 -*-
import os
import numpy as np
from numpy import sqrt
import scipy.integrate as integrate
import sys
import csv
import pandas as pd
import openpyxl
from matplotlib import pyplot as plt
from matplotlib import gridspec
from tklib.tkutils import terminate, pfloat, pint, getarg, getintarg, getfloatarg
from tklib.tkapplication import tkApplication
from tklib.tkvariousdata import tkVariousData
from tklib.tktransport.tkTransport import LorentzNumber_FEA, FermiIntegral_x2
from tklib.tktransport.tkTransport import FermiIntegral_fast as FermiIntegral
from tklib.tktransport.tkTransport import meff2NC_FEA, NC2meff_FEA, meff2DC0_FEA
#===================================
# physical constants
#===================================
pi = 3.14159265358979323846
e = 1.60218e-19 # C";
kB = 1.380658e-23 # JK<sup>-1</sup>";
kBe = kB / e #8.6173303e-5 # kB/e
kpi05 = 2.0 / sqrt(pi)
# Global variables
infile = "20221219Ba3GeO.xlsx"
outfile = None
meff = 1.0
T_label = 0
S_label = 1
sigma_label = 2
kappa_label = 3
nmaxiter = 100
eps = 1.0e-10
a = 1000 #max
b = -300 #min
#===================================
# figure configuration
#===================================
figsize = (12, 8)
fontsize = 16
legend_fontsize = 12
markersize = 8.0
app = None
#===================================
# Treat arguments
#===================================
infile = getarg(1, defval = infile)
T_label = getarg(2, defval = T_label)
S_label = getarg(3, defval = S_label)
sigma_label = getarg(4, defval = sigma_label)
kappa_label = getarg(5, defval = kappa_label)
meff = getfloatarg(6, defval = meff)
header, ext = os.path.splitext(infile)
filebody = os.path.basename(header)
outfile = f'{header}-out.xlsx'
#===================================
# Other functions
#===================================
[ドキュメント]
def usage(app):
"""
スクリプトの正しい使用方法を表示します。
実行時の引数の指定方法に関するガイダンスを標準出力に表示します。
:param app: アプリケーションインスタンス (tkApplication)。エラー終了時に使用されることがあります。
:type app: tklib.tkapplication.tkApplication
:returns: なし
:rtype: None
"""
argv = sys.argv
print("")
print("Usage: python {} infile T_label S_label(S/cm) sigma_label(S/cm) kappa_label(W/m/K)".format(argv[0]))
print(" ex: python {} single '{}' '{}' '{}' '{}' '{}'".format(argv[0],
infile, T_label, S_label, sigma_label, kappa_label))
[ドキュメント]
def savecsv(outfile, header, datalist):
"""
データをCSVファイルに保存します。
指定されたファイル名でCSVファイルを新規作成し、ヘッダーとデータリストを書き込みます。
書き込み中にエラーが発生した場合は、エラーメッセージを標準出力に表示します。
:param outfile: 出力CSVファイルのパス。
:type outfile: str
:param header: CSVのヘッダー行のリスト。
:type header: list[str]
:param datalist: 保存するデータを含むリストのリスト。各サブリストは列を表します。
:type datalist: list[list[float]]
:returns: なし
:rtype: None
"""
try:
print("Write to [{}]".format(outfile))
f = open(outfile, 'w')
except:
# except IOError:
print("line 98 Error: Can not write to [{}]".format(outfile))
else:
fout = csv.writer(f, lineterminator='\n')
fout.writerow(header)
# fout.writerows(data)
for i in range(0, len(datalist[0])):
a = []
for j in range(len(datalist)):
a.append(datalist[j][i])
fout.writerow(a)
f.close()
[ドキュメント]
def read_file(fname):
"""
指定されたファイルからデータを読み込みます。
tkVariousDataクラスを使用して、指定されたファイルからラベルとデータリストを抽出します。
ファイル形式はtkVariousDataがサポートするものに従います。
:param fname: 読み込むファイルのパス。
:type fname: str
:returns: (データファイルオブジェクト, ラベルリスト, データリスト) のタプル。
:rtype: tuple[tklib.tkvariousdata.tkVariousData, list[str], list[list[float]]]
"""
print("")
datafile = tkVariousData(infile)
labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = usage)
return datafile, labels, datalist
#Calculation method
[ドキュメント]
def cal_ke(sigma, L, T):
"""
電子熱伝導率 (kappa_e) を計算します。
Wiedemann-Franzの法則に基づいて、電気伝導率、ローレンツ数、および温度から電子熱伝導率を計算します。
この関数では、入力される電気伝導率sigmaの単位はS/mを想定しています。
:param sigma: 電気伝導率 (S/m)。
:type sigma: float
:param L: ローレンツ数 (WΩ/K²)。
:type L: float
:param T: 温度 (K)。
:type T: float
:returns: 電子熱伝導率 (W/m/K)。
:rtype: float
"""
return sigma * L * T
[ドキュメント]
def Fr(EF, r):
"""
フェルミ・ディラック積分 (F_r(eta)) を計算します。
積分下限0から無限大までのE^r / (1 + exp(E - EF))の積分を計算します。
ここでEは規格化されたエネルギー、EFは規格化されたフェルミ準位 (eta = E_F/k_B T) です。
`scipy.integrate.quad`を使用して数値積分を実行します。
:param EF: 規格化されたフェルミ準位 (eta = E_F/k_B T)。
:type EF: float
:param r: フェルミ積分の次数。
:type r: float
:returns: 計算されたフェルミ・ディラック積分値。
:rtype: float
"""
result = integrate.quad(lambda E: E**r / (1 + np.exp(E - EF)), 0, np.inf)
#print("felmi_integral = {}".format(result))
fel = float(result[0])
#print(fel)
return fel
[ドキュメント]
def S_comparison(T, S, sigma, EF):
"""
測定されたゼーベック係数とSPBモデルから計算されたゼーベック係数の差を計算します。
Single Parabolic Band (SPB) モデルに基づき、与えられた規格化フェルミ準位 (EF) で計算されるゼーベック係数と、
実験的に測定されたゼーベック係数 (S) の差を評価します。この関数は通常、
フェルミ準位を決定するための数値解法（例：二分法）内で使用されます。
フェルミ・ディラック積分は、次数 r=0.0 および r=1.0 で計算されます。
引数`sigma`は渡されますが、この関数内では直接使用されません。
:param T: 温度 (K)。
:type T: float
:param S: 実験的に測定されたゼーベック係数 (V/K)。
:type S: float
:param sigma: 電気伝導率 (S/cmまたはS/m)。この関数では直接使用されません。
:type sigma: float
:param EF: 規格化されたフェルミ準位 (eta = E_F/k_B T)。
:type EF: float
:returns: 測定SとSPBモデルから計算されたSの差 (V/K)。
:rtype: float
"""
diff = S - kBe * (2.0 * Fr(EF, r = 1.0) / Fr(EF, r = 0.0) - EF)
#print("diff = {}".format(diff))
return diff
[ドキュメント]
def Lorenz(EF):
"""
規格化フェルミ準位からローレンツ数を計算します。
Single Parabolic Band (SPB) モデルに基づき、規格化されたフェルミ準位 (EF/kT) から
ローレンツ数 (L) を計算します。フェルミ・ディラック積分は、次数 r=0.0, r=1.0, r=2.0 で使用されます。
:param EF: 規格化されたフェルミ準位 (eta = E_F/k_B T)。
:type EF: float
:returns: 計算されたローレンツ数 (WΩ/K²)。
:rtype: float
"""
answer = kBe**2 * (3.0 * Fr(EF, r = 2.0) / Fr(EF, r = 0.0) - (2.0 * Fr(EF, r = 1.0) / Fr(EF, r = 0.0))**2)
return answer
[ドキュメント]
def cal_EF_from_S_bisection(T, S, sigma, max_ef, min_ef, err):
"""
二分法を用いてゼーベック係数から規格化フェルミ準位 (EF/kT) を推定します。
与えられた温度 (T) とゼーベック係数 (S) に対して、`S_comparison`関数を評価しながら
指定された探索範囲 (min_ef, max_ef) で二分法を実行し、規格化フェルミ準位 (EF/kT) を収束させます。
収束しない場合、エラーメッセージを出力し、プログラムを終了します。
引数`sigma`は`S_comparison`関数に渡されますが、この関数内では直接使用されません。
この関数に渡される`S`は絶対値が取られ、`sigma`はS/mに変換されて渡されることが想定されます。
:param T: 温度 (K)。
:type T: float
:param S: 測定されたゼーベック係数 (V/K)。この関数に渡される前に絶対値が取られます。
:type S: float
:param sigma: 電気伝導率 (S/m)。`S_comparison`関数に渡されますが、この関数では直接使用されません。
:type sigma: float
:param max_ef: フェルミ準位の探索範囲の上限（`a`に相当する初期値、例: 1000）。
:type max_ef: float
:param min_ef: フェルミ準位の探索範囲の下限（`b`に相当する初期値、例: -300）。
:type min_ef: float
:param err: 収束判定に使用される許容誤差。
:type err: float
:returns: 推定された規格化フェルミ準位 (EF/kT)。
:rtype: float
:raises SystemExit: 探索範囲に解がない場合、または最大反復回数を超えて収束しなかった場合。
"""
xEF = (max_ef + min_ef) / 2.0
ymin = S_comparison(T, S, sigma, min_ef)
ymax = S_comparison(T, S, sigma, max_ef)
yh = S_comparison(T, S, sigma, xEF)
print(f"xEF={xEF} ymin={ymin} ymax={ymax} yh={yh} min={min_ef} max={max_ef}")
if ymin * ymax > 0.0:
print("line Error: ymin * ymax > 0.0")
app.terminate(pause = True)
for i in range(nmaxiter):
print(f"xEF={xEF:8.3f} d=({ymin:8.4g} {yh:8.4g} {ymax:8.4g}) min={min_ef:8.3f} max={max_ef:8.3f}")
if abs(yh) < eps:
# if abs(yh) < eps or abs(max_ef - min_ef) < eps:
return xEF
if ymin * yh > 0.0:
min_ef = xEF
else:
max_ef = xEF
xEF = (max_ef + min_ef) / 2.0
ymin = S_comparison(T, S, sigma, min_ef)
ymax = S_comparison(T, S, sigma, max_ef)
yh = S_comparison(T, S, sigma, xEF)
else:
print(f"Did not converge in {nmaxiter} iterations.")
exit()
[ドキュメント]
def cal_Ne(T, xEF, NC):
"""
キャリア濃度 (Ne) を計算します。
Single Parabolic Band (SPB) モデルに基づき、温度 (T)、規格化フェルミ準位 (xEF = EF/kT)、
および有効状態密度 (NC) からキャリア濃度を計算します。
フェルミ・ディラック積分は次数 r=0.5 で使用されます。
引数`T`は渡されますが、この関数内では直接使用されません。
:param T: 温度 (K)。この関数では直接使用されません。
:type T: float
:param xEF: 規格化されたフェルミ準位 (eta = E_F/k_B T)。
:type xEF: float
:param NC: 有効状態密度 (cm^-3)。tklib.tktransport.tkTransport.meff2NC_FEAで計算されます。
:type NC: float
:returns: 計算されたキャリア濃度 (cm^-3)。
:rtype: float
"""
F05 = Fr(xEF, r = 0.5)
N = NC * kpi05 * F05
return N
[ドキュメント]
def main():
"""
スクリプトのメイン実行関数です。
入力ファイルの読み込み、ゼーベック係数からのフェルミ準位、ローレンツ数、
電子熱伝導率の計算、結果のファイルへの保存、およびグラフの描画を行います。
コマンドライン引数から入力ファイル名とデータ列のラベルを読み取ります。
"""
global app
app = tkApplication()
logfile = app.replace_path(infile)
print(f"Open logfile [{logfile}]")
app.redirect(targets = ["stdout", logfile], mode = 'w')
print("")
print("Estimate Lorentz factor from Seebeck coefficient with r = 0.0")
print("input file :", infile)
print("output file:", outfile)
print("T_label :", T_label)
print("S_label :", S_label)
print("sigma_label:", sigma_label)
print("kappa_label:", kappa_label)
print("")
if '***' in T_label:
app.terminate("Error: Choose T", usage = usage, pause = True)
if '***' in S_label:
app.terminate("Error: Choose S", usage = usage, pause = True)
if '***' in sigma_label:
app.terminate("Error: Choose sigma", usage = usage, pause = True)
if '***' in kappa_label:
app.terminate("Error: Choose kappa", usage = usage, pause = True)
if 'S/m' in sigma_label or 'cm' not in sigma_label:
app.terminate(f"Error: The unit of sigma must be 'S/cm'. The given label is [{sigma_label}]", usage = usage, pause = True)
if 'uK' in S_label or 'microK' in S_label:
app.terminate(f"Error: The unit of S must be 'V/K'. The given label is [{S_label}]", usage = usage, pause = True)
print("Read data from [{}]".format(infile))
datafile, labels, datalist = read_file(infile)
T_idx = datafile.FindLabelIndex(T_label, flag = 'i')
S_idx = datafile.FindLabelIndex(S_label, flag = 'i')
sigma_idx = datafile.FindLabelIndex(sigma_label, flag = 'i')
kappa_idx = datafile.FindLabelIndex(kappa_label, flag = 'i')
print("T index: ", T_idx)
print("S index: ", S_idx)
print("sigma index: ", sigma_idx)
print("kappa index: ", kappa_idx)
Tlabel, Tlist = datafile.FindDataArray(T_label, flag = 'i')
Slabel, Slist = datafile.FindDataArray(S_label, flag = 'i') # V/K
sigmalabel, sigmalist = datafile.FindDataArray(sigma_label, flag = 'i') # S/cm
kappalabel, kappalist = datafile.FindDataArray(kappa_label, flag = 'i')
print("")
print(f"meff = {meff} me")
print(f"{'T(K)':8} {'S(V/K)':10} {'sigma(S/cm)':10} {'kappa':10}")
for i in range(len(Tlist)):
print(f"{Tlist[i]:8.3g} {Slist[i]:10.4g} {sigmalist[i]:10.4g} {kappalist[i]:10.4g}")
print("")
print(f"meff = {meff} me")
EFlist = []
NClist = []
Nelist = []
EFkTlist = []
Llist = []
kelist = []
for ic in range(len(Tlist)):
T = Tlist[ic]
S = Slist[ic]
sigma = sigmalist[ic]
kappa = kappalist[ic]
print(f"{ic:3d}: T={T} K S={S} V/K sigma={sigma} S/cm kappa={kappa}")
# Sの絶対値を使用し、sigmaをS/mに変換して渡す
xEF = cal_EF_from_S_bisection(T, abs(S), sigma * 1.0e-2, max_ef = a, min_ef = b, err = 1e-10)
EF = xEF * kB * T / e
NC = meff2NC_FEA(meff, T) # cm-3
Ne = cal_Ne(T, xEF, NC) # cm-3
L = Lorenz(xEF)
# sigmaをS/mに変換してcal_keに渡す
ke = cal_ke(sigma * 1.0e-2, L, T) # sigma to S/m
print(f"T = {T} K")
print(f"S = {S*1.0e6} uV/K")
print(f"EF/kT = {xEF}")
print(f"EF = {EF} eV")
print(f"NC = {NC} cm^-3")
print(f"Ne = {Ne} cm^-3")
print(f"L = {L}")
print(f"ke = {ke} W/m/K")
print(f"kT/e = { kB * T / e}")
EFlist.append(EF)
EFkTlist.append(xEF)
NClist.append(NC)
Nelist.append(Ne)
Llist.append(L)
kelist.append(ke)
df = pd.DataFrame(np.array([Tlist, Nelist, Slist, sigmalist, kappalist, EFkTlist, EFlist, Llist, kelist]).T,
columns = ["T(K)", "Ne(cm^-3)", "S(V/K)", "simga(S/cm)", "kappa(W/m/K)", "EF/kBT", "EF(eV)", "L", "kappa_e(W/m/K)"])
print("")
print("Save L data to [{}]".format(outfile))
df.to_excel(outfile, index = False, header = True)
#=============================
# グラフの表示
#=============================
print("")
fig = plt.figure(figsize = figsize)
ax1 = fig.add_subplot(2, 2, 1)
ax1b = ax1.twinx()
ax2 = fig.add_subplot(2, 2, 3)
ax3 = fig.add_subplot(2, 2, 2)
ax3b = ax3.twinx()
ax4 = fig.add_subplot(2, 2, 4)
ax1.tick_params(labelsize = fontsize)
ax2.tick_params(labelsize = fontsize)
ax3.tick_params(labelsize = fontsize)
ax4.tick_params(labelsize = fontsize)
ins1 = ax1.plot(Tlist, Nelist, label = '$N_e$ (cm$^{-3}$)', linestyle = '', marker = 'o', markerfacecolor = 'black', markeredgecolor = 'black')
ax1.set_xlabel("$T$ (K)", fontsize = fontsize)
ax1.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
ax1.set_yscale('log')
S = [Slist[i] * 1.0e6 for i in range(len(Slist))]
ins2 = ax1b.plot(Tlist, S, label = '$S$ ($\\mu$V/K)', linestyle = '', marker = 's', markerfacecolor = 'blue', markeredgecolor = 'blue')
# ax1b.set_xlabel("$T$ (K)", fontsize = fontsize)
ax1b.set_ylabel("$S$ ($\\mu$V/K)", fontsize = fontsize)
ins = ins1 + ins2
ax1.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
ax2.plot(Tlist, Llist, label = '$L$', linestyle = '', marker = 's')
ax2.set_xlabel("$T$ (K)", fontsize = fontsize)
ax2.set_ylabel("$L$", fontsize = fontsize)
ins1 = ax3.plot(Nelist, sigmalist, label = '$\\sigma$ (S/cm)', linestyle = '', marker = 'o', markerfacecolor = 'black', markeredgecolor = 'black')
ax3.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
ax3.set_ylabel("$\\sigma$ (S/cm)", fontsize = fontsize)
ax3.set_xscale('log')
ins2 = ax3b.plot(Nelist, S, label = '$S$ ($\\mu$V/K)', linestyle = '', marker = 's', markerfacecolor = 'blue', markeredgecolor = 'blue')
ax3b.set_ylabel("$S$ ($\\mu$V/K)", fontsize = fontsize)
ins = ins1 + ins2
ax3.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
ax4.plot(Nelist, Llist, label = '$L$', linestyle = '', marker = 's')
ax4.set_xlabel("$N_e$ (K)", fontsize = fontsize)
ax4.set_ylabel("$L$", fontsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
print()
print("NOTE: r is limited to 0")
print()
app.terminate("", usage = usage, pause = True)
if __name__ == '__main__':
main()