"""
Seebeck2.py
============
概要:
ゼーベック係数と温度から有効質量を推定するためのスクリプト。
熱電材料の電気伝導度、ゼーベック係数、ホールキャリア濃度、ホール移動度などの
輸送特性を計算、シミュレーション、およびフィッティングを行う。
詳細説明:
このスクリプトは、G.J. Snyder, A. Pereyra and R. Gurunathan, Adv. Funct. Mater. 32, 2112772 (2022)
に記載されている有効質量推定手法を実装している。
半導体の輸送特性を計算するために、Fermi-Dirac統計とさまざまな散乱メカニズムを考慮したモデルを使用する。
以下のモードをサポートしている:
- basic: フェルミ・ディラック分布とフェルミ積分の基本関数をプロットする。
- prop: ゼーベック係数、電気伝導度、移動度、ローレンツ数、ZT値などの熱電特性をプロットする。
- Hall: ホール係数、ホール因子、ホール移動度などのホール関連特性をプロットする。
- init: 最適化のための初期パラメータファイル（.inファイル）を生成する。
- sim: 入力ファイルからデータを読み込み、熱電特性をシミュレーションし、ヨンカープロットやピサレンコプロットを表示する。
- fit: 実測のゼーベック係数データと電気伝導度データにフィッティングを行う。
- T: 温度依存性を考慮した熱電特性のシミュレーションとプロットを行う。
- calL: 実測の電気伝導度、キャリア濃度、温度からローレンツ数と電子熱伝導率を計算する。
- calF: 実測のホールキャリア濃度とホール移動度からホール因子を計算する。
関連リンク:
- `G.J. Snyder, A. Pereyra and R. Gurunathan, Adv. Funct. Mater. 32, 2112772 (2022)
<https://onlinelibrary.wiley.com/doi/10.1002/adfm.202112772>`_
- :doc:`Seebeck2_usage`
"""
import os
import sys
import openpyxl
from math import exp, sqrt, log, gamma
import numpy as np
#from numpy import arange
from scipy import optimize # newton関数はscipy.optimizeモジュールに入っている
from scipy.interpolate import interp1d
from matplotlib import pyplot as plt
from scipy import optimize # newton関数はscipy.optimizeモジュールに入っている
from scipy.optimize import minimize
from tklib.tkutils import terminate, getarg, getintarg, getfloatarg, safe_getelement, validate_error
from tklib.tkutils import pint, pfloat, sort_lists
from tklib.tkinifile import tkIniFile
from tklib.tkvariousdata import tkVariousData
from tklib.tksci.tksci import pi, h, hbar,c, e, kB, me
from tklib.tksci.tkmatrix import make_matrix1, make_matrix2
from tklib.tkapplication import tkApplication
from tklib.tktransport.tkTransport import read_datafile, FermiIntegral_x2
from tklib.tktransport.tkTransport import FermiIntegral_fast as FermiIntegral
from tklib.tktransport.tkTransport import integrate_Simpson_list, fe, fh, meff2NC_FEA, meff2DC0_FEA, BMShift_FEA
from tklib.tktransport.tkWeightedMobility import weighted_mobility, weighted_mobility_new, weighted_mobility_exact
from tklib.tktransport.tkTransport import LorentzNumber_FEA
from tklib.tktransport.tkDOS_FEA import integrate_Simpson, integrate_Simpson_list, tkDOS
from tklib.tktransport.tkmobility_tau import tkMobility, split_optstr
"""
Estimate m* from S and T
# G.J. Snyder, A. Pereyra and R. Gurunathan, Adv. Funct. Mater. 32, 2112772 (2022)
# Effective mass from Seebeck coefficient
# https://onlinelibrary.wiley.com/doi/10.1002/adfm.202112772
"""
#================================
# Global variables
#================================
Debug = 0
L_FEA = LorentzNumber_FEA
#mode: 'basic', 'prop', 'Hall', 'init', 'sim', 'fit', 'T', 'calL', 'calF'
mode = 'sim'
infile = 'SnSeTe-S-Hall.xlsx'
T_label = r'^T[\s|\(|\[|$]'
sigma_label = r'^sigma[$|\S]?.*$'
n_label = r'^N[$|\S]?.*$'
mu_label = r'^mu[$|\s|\(\[]?.*$'
S_label = r'^S[$|\s|\(\[]?.*$'
outxlsFjfile = "Fj.xlsx"
outxlsfFDfile = "fFD.xlsx"
outxlsPropfile = "properties.xlsx"
outxlsSLfile = "S-L.xlsx"
outcsvFitfile = ""
outcsvLfile = ""
parameterfile = None
parameterbkfile = None
outxlsfile = None
# Calculation ranges for mode = 'basic'
# x = E / kB / T
xmin_basic = -20.0
xmax_basic = 20.0
nx_basic = 401
xstep_basic = (xmax_basic - xmin_basic) / (nx_basic - 1)
# Calculation ranges for mode = 'prop'
# x = E / kB / T
xmin = -30.0
xmax = 250.0
nx = 281
xstep = (xmax - xmin) / (nx - 1)
# Seebeck coefficient
Smin = 0.0
Smax = 1.0e-3 # V/K
nS = 101
Sstep = (Smax - Smin) / (nS - 1)
# N range
Nmin = 1.0e15 # cm-3
Nmax = 1.0e22 # cm-3
nN = 101
# T range
T0 = 300.0 # K
Tmin = 300.0
Tmax = 800.0
nT = 6
# read from DOSCAR
dos = tkDOS()
dos.meeff = 0.3
#dos.mheff = 1.0
#dos.EV = 0.0
dos.EC = 0.0
#dos.EA = 0.05 # Acceptor level, eV
#dos.NA = 0.0e17
#dos.ED = 1.05 # Donor level
#dos.ND = 0.0e17
#dos.EF0 = 0.0 # initial EF to find EF
#移動度パラメータ
mobility = tkMobility()
# For debug purpose. 1 will use 3-terms polynomial for the inverse of mobility
mobility.debug = 0
mobility.use_simple = 0
#mobility.use_simple = 1
# charge, effective mass, scattering factor
mobility.charge = 1.0 # in e, 1.0 for hole, -1.0 for electron
mobility.meff = dos.meeff # in me
mobility.rfac = 0.5 # tau = (meff/2)^0.5 * l0(T) * E^(r-0.5)
mobility.l0 = 1.0e-8 # m
# Lattice thermal conductivity
klatt = 5.0 # W/m/K
# EF range for mode = 'EF'
#dEFmin = -1.0 # measured from EV, eV
#dEFmax = 1.0 # measured from EC, eV
#dos.nEF = 50
#dos.Estep = 0.01 # Integration step, eV
epsEF = 1.0e-5
label_sample = None
xsample = None
label_S = None
yS = None
label_sigma = None
ysigma = None
label_N = None
yN = None
label_mu = None
ymu = None
ysigmaini = None
ymuini = None
ySini = None
# 最適化パラメータの初期値
varname = ["meff", "l0"]
varunit = [ "me", "m"]
ai0 = []
optid = [ 1, 1]
app = None
#=============================================
# scipy.optimize.minimizeで使うアルゴリズム
#=============================================
#nelder-mead Downhill simplex
#powell Modified Powell
#cg conjugate gradient (Polak-Ribiere method)
#bfgs BFGS法
#newton-cg Newton-CG
#trust-ncg 信頼領域 Newton-CG 法
#dogleg 信頼領域 dog-leg 法
#L-BFGS-B’ (see here)
#TNC’ (see here)
#COBYLA’ (see here)
#SLSQP’ (see here)
#trust-constr’(see here)
#dogleg’ (see here)
#trust-exact’ (see here)
#trust-krylov’ (see here)
method = "nelder-mead"
#method = 'cg'
#method = 'powell'
#method = 'bfgs'
maxiter = 1000
tol = 1.0e-4
h_diff = 1.0e-3
outputinterval = 1
Nmin_fit = '*'
Nmax_fit = '*'
sigmamin_fit = '*'
sigmamax_fit = '*'
#=============================
# Graph configuration
#=============================
fig = None
figsize = (12, 8)
figsize_sim = (10, 6)
figsize_small = (8, 6)
figsize_FD = (8, 6)
fontsize = 18
legend_fontsize = 12
graphupdateinterval = 10
#=============================
# Treat argments
#=============================
[ドキュメント]
def parameter_list():
"""
概要:
現在の移動度パラメータのリストを返す。
戻り値:
tuple: (mobility.meff, mobility.l0) のタプル。
"""
return mobility.meff, mobility.l0
[ドキュメント]
def set_parameters(ai):
"""
概要:
与えられた値で移動度とDOSのパラメータを更新する。
詳細説明:
`ai` リストの最初の要素を有効質量 `dos.meeff` および `mobility.meff` に設定する。
もし `ai` に2つ以上の要素があれば、2番目の要素を平均自由行程のプレファクター `mobility.l0` に設定し、
その後、散乱パラメータを更新する。
引数:
:param ai: (list) 設定するパラメータ値のリスト。
`ai[0]` は有効質量 `meff` (単位 `me`)。
`ai[1]` は平均自由行程のプレファクター `l0` (単位 `m`)。
"""
global dos, mobility
dos.meeff = ai[0]
dos.NC = meff2NC_FEA(dos.meeff, T0)
dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
mobility.meff = ai[0]
if len(ai) > 1:
# print("ai1=", ai[1])
mobility.l0 = ai[1]
mobility.set_scattering_parameters()
# print("l0=", mobility.l0)
[ドキュメント]
def read_parameters(path):
"""
概要:
指定されたパスのINIファイルからパラメータを読み込む。
詳細説明:
INIファイルから `meff`, `l0`, `charge`, `rfac`, `EC`, `T0` などの設定を読み込み、
対応するグローバル変数（`ai0`, `optid`, `mobility`, `dos`, `T0`）を更新する。
読み込みに失敗した場合は何もせずに終了する。
引数:
:param path: (str) パラメータが記述されたINIファイルのパス。
"""
global ai0, optid, T0
ini = tkIniFile()
inf = ini.ReadAll(path, AddSection = 0)
if inf is None:
return
ai0 = list(ai0)
keylist = ["meff", "l0"]
for i in range(2):
key = keylist[i]
str = inf.get(key, None)
val, id = split_optstr(str)
if val is not None:
ai0[i] = val
optid[i] = id
mobility.meff = pfloat(ai0[0])
dos.meeff = pfloat(ai0[0])
mobility.l0 = pfloat(ai0[1])
mobility.charge = pfloat(safe_getelement(inf, "charge", mobility.charge))
mobility.rfac = pfloat(safe_getelement(inf, "rfac", mobility.rfac))
dos.EC = pfloat(safe_getelement(inf, "EC", dos.EC))
T0 = pfloat(safe_getelement(inf, "T0", T0))
[ドキュメント]
def save_parameters(path, ai, args):
"""
概要:
現在のパラメータをINIファイルに保存する。
詳細説明:
`Preferences` セクションには `args` 辞書の内容を、`Parameters` セクションには
最適化パラメータ `ai` とそれに対応する `optid` を保存する。
引数:
:param path: (str) パラメータを保存するINIファイルのパス。
:param ai: (list) 最適化パラメータのリスト。
:param args: (dict) その他の設定を格納した辞書。
"""
ini = tkIniFile(path)
for key in args.keys():
ini.WriteString('Preferences', key, args[key])
for i in range(len(ai)):
ini.WriteString('Parameters', varname[i], "{}:{}".format(ai[i], optid[i]))
[ドキュメント]
def save_parameterfile(ai = ai0, S2 = ''):
"""
概要:
現在のパラメータと計算結果をINIファイルに保存する。
詳細説明:
グローバル変数 `parameterfile` を使用して、入力ファイル名、基準温度、
キャリア電荷、散乱因子、伝導帯端エネルギー、有効状態密度、状態密度定数、
およびフィット結果の二乗誤差 `S2` をINIファイルに書き込む。
引数:
:param ai: (list, optional) 保存するパラメータのリスト。デフォルトはグローバル変数 `ai0`。
:param S2: (str, optional) フィットの二乗誤差。デフォルトは空文字列。
"""
save_parameters(parameterfile, ai,
{"infile": infile, "T0": T0,
"charge": mobility.charge, "rfac": mobility.rfac,
"EC": dos.EC, "NC": dos.NC, "DC0": dos.DC0,
"S2": S2})
[ドキュメント]
def print_parameters(ai = None):
"""
概要:
現在のパラメータをコンソールに出力する。
詳細説明:
グローバル変数 `varname`, `varunit`, `optid` を使用して、
各パラメータの名前、値、単位、および最適化IDを表示する。
引数:
:param ai: (list, optional) 表示するパラメータのリスト。デフォルトはグローバル変数 `ai0`。
"""
if ai is None:
ai = ai0
for i in range(len(varname)):
print(" {:10s}: {:14.8g} {:6} optid={}".format(varname[i], ai[i], varunit[i], optid[i]))
[ドキュメント]
def usage(app):
"""
概要:
スクリプトの利用方法とコマンドライン引数の例を表示する。
詳細説明:
各実行モード（`basic`, `prop`, `Hall`, `init`, `sim`, `fit`, `T`, `calL`, `calF`）の
目的と必要なコマンドライン引数のフォーマットを説明する。
引数:
:param app: (tkApplication) アプリケーションオブジェクト。
"""
print("")
print("Usage: Variables in () are optional")
print(" (i) python {} basic".format(sys.argv[0]))
print(" Plot basic functions")
print(" (ii) python {} prop T meff r l0 kappa_latt".format(sys.argv[0]))
print(" meff in me0, l0 in m, kappa_latt in W/m/K")
print(" Plot Seebeck related properties")
print(" ex: python {} basic {} {} {} {} {}".format(sys.argv[0], T0, dos.meeff, mobility.rfac, mobility.l0, klatt))
print(" (ii') python {} Hall T meff r l0".format(sys.argv[0]))
print(" Plot Hall related properties")
print(" ex: python {} Hall {} {} {} {}"
.format(sys.argv[0], T0, dos.meeff, mobility.rfac, mobility.l0))
print(" (iii) python {} init infile".format(sys.argv[0]))
print(" Create .in file")
print(" (iv) python {} sim infile T meff r l0".format(sys.argv[0]))
print(" Calculate weighted mobility etc from input file")
print(" Plot Jonker / Pisarenko plots")
print(" ex: python {} sim {} {} {} {} {}".format(sys.argv[0], infile, T0, dos.meeff, mobility.rfac, mobility.l0))
print(" (v) python {} fit T meff r l0".format(sys.argv[0]))
print(" Fit to Jonker / Pisarenko plots")
print(" ex: python {} fit {} {} {} {} {}"
.format(sys.argv[0], infile, T0, mobility.meff, mobility.rfac, mobility.l0))
print(" (vi) python {} T infile Tmin Tmax nT".format(sys.argv[0]))
print(" Simulate T dependences")
print(" ex: python {} fit {} {} {} {}"
.format(sys.argv[0], infile, Tmin, Tmax, nT))
print(" (vii) python {} calL infile meff r l0 ".format(sys.argv[0]))
print(" Calculate L and kappa,e from T, Ne and sigma")
print(" ex: python {} calL {} {} {} {}"
.format(sys.argv[0], infile, dos.meeff, mobility.rfac, mobility.l0))
print(" (viii) python {} calF infile meff r l0 T_label n_label mu_label sigma_label".format(sys.argv[0]))
print(" Calculate FHall from T, Ne and mu")
print(" ex: python {} calF {} {} {} {} {} {} {}"
.format(sys.argv[0], infile, dos.meeff, mobility.rfac, mobility.l0,
T_label, n_label, mu_label, sigma_label))
[ドキュメント]
def updatevars():
"""
概要:
コマンドライン引数を解析し、グローバル変数を更新する。
詳細説明:
`sys.argv` を参照し、実行モードに応じて以下のグローバル変数を更新する:
`mode`, `infile`, `outxlsxfile`, `parameterfile`, `parameterbkfile`,
`T_label`, `n_label`, `mu_label`, `sigma_label`, `S_label`,
`mobility`, `dos`, `klatt`, `T0`, `Tmin`, `Tmax`, `nT`,
`ai0`, `optid`, `method`, `tol`, `maxiter`, `h_diff`,
`Nmin_fit`, `Nmax_fit`, `sigmamin_fit`, `sigmamax_fit`。
不正なモードが指定された場合は、`usage` メッセージを表示して終了する。
引数:
なし。
"""
global mode, infile, outxlsxfile, parameterfile, parameterbkfile
global T_label, n_label, mu_label, sigma_label, S_label
global mobility, dos, klatt
global T0, Tmin, Tmax, nT
global ai0, optid
global method, tol, maxiter, h_diff
global Nmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fit
argv = sys.argv
# if len(argv) == 1:
# usage()
# exit()
mode = getarg( 1, mode)
if mode == 'basic':
pass
elif mode == 'prop':
T0 = getfloatarg(2, T0)
dos.meeff = getfloatarg(3, dos.meeff)
mobility.rfac = getfloatarg(4, mobility.rfac)
mobility.l0 = getfloatarg(5, mobility.l0)
klatt = getfloatarg(6, klatt)
elif mode == 'Hall':
T0 = getfloatarg(2, T0)
dos.meeff = getfloatarg(3, dos.meeff)
mobility.rfac = getfloatarg(4, mobility.rfac)
mobility.l0 = getfloatarg(5, mobility.l0)
elif mode == 'init':
infile = getarg (2, infile)
elif mode == 'sim':
infile = getarg ( 2, infile)
T_label = getarg ( 3, T_label)
n_label = getarg ( 4, n_label)
mu_label = getarg ( 5, mu_label)
sigma_label = getarg ( 6, sigma_label)
S_label = getarg ( 7, S_label)
T0 = getfloatarg( 8, T0)
dos.meeff = getfloatarg( 9, dos.meeff)
mobility.rfac = getfloatarg(10, mobility.rfac)
mobility.l0 = getfloatarg(11, mobility.l0)
elif mode == 'fit':
infile = getarg (2, infile)
T_label = getarg (3, T_label)
n_label = getarg (4, n_label)
mu_label = getarg (5, mu_label)
sigma_label = getarg (6, sigma_label)
S_label = getarg (7, S_label)
T0 = getfloatarg(8, T0)
dos.meeff = getfloatarg(9, dos.meeff)
mobility.rfac = getfloatarg(10, mobility.rfac)
mobility.l0 = getfloatarg(11, mobility.l0)
method = getarg (12, method)
tol = getfloatarg(13, tol)
maxiter = getintarg (14, maxiter)
Nmin_fit = getarg (15, Nmin_fit)
Nmax_fit = getarg (16, Nmax_fit)
sigmamin_fit = getarg (17, sigmamin_fit)
sigmamax_fit = getarg (18, sigmamax_fit)
elif mode == 'T':
infile = getarg (2, infile)
Tmin = getfloatarg(3, Tmin)
Tmax = getfloatarg(4, Tmax)
nT = getintarg (5, nT)
elif mode == 'calL' or mode == 'calF':
infile = getarg (2, infile)
T_label = getarg (3, T_label)
n_label = getarg (4, n_label)
mu_label = getarg (5, mu_label)
sigma_label = getarg (6, sigma_label)
dos.meeff = getfloatarg(7, dos.meeff)
mobility.rfac = getfloatarg(8, mobility.rfac)
mobility.l0 = getfloatarg(9, mobility.l0)
elif mode == 'help' or mode == 'usage':
app.terminate("", usage = usage, pause = True)
else:
app.terminate("Error in updatevars(): Invalid mode {}".format(mode), usage = usage, pause = True)
mobility.meff = dos.meeff
mobility.set_scattering_parameters(mobility.l0, mobility.rfac)
ai0 = parameter_list()
[ドキュメント]
def basic():
"""
概要:
フェルミ・ディラック分布とフェルミ積分に関する基本的なプロットとデータ出力を行う。
詳細説明:
`xmin_basic` から `xmax_basic` の範囲で、フェルミ・ディラック分布関数 `fFD`、
そのホール版 `fFDh`、微分値 `-df/dx`、および近似式 `fFDa` を計算する。
また、フェルミ積分 `F_r` を様々な `r` (0, 0.5, 1, 1.5, 2) について計算し、
結果をExcelファイル (`outxlsfFDfile`, `outxlsFjfile`) に出力し、
matplotlibでグラフ表示する。
引数:
なし。
"""
global mobility, dos
print("")
print("mode:", mode)
print("")
print("Fermi-Dirac distribution")
xx = []
yfFD = []
yfFDh = []
ymdfdx = []
yfFDapprox = []
print("{:8}\t{:12}\t{:12}\t{:12}\t{:12}".format("x=E/kBT", "fFDe(exact)", "fFDe(approx)", "fFDh(exact)", "-df/dx"))
for i in range(nx_basic):
x = xmin_basic + i * xstep_basic
xx.append(x)
fFD = 1.0 / (exp(x) + 1.0)
fFDh = 1.0 / (exp(-x) + 1.0)
mdfdx = fFD * fFDh
fFDa = 0.5 - 0.25 * x + 1.0 / 48.0 * x**3 - 1.0 / 64.0 * x**5
yfFD.append(fFD)
yfFDh.append(fFDh)
ymdfdx.append(mdfdx)
yfFDapprox.append(fFDa)
print("{:8.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}".format(x, fFD, fFDa, fFDh, mdfdx))
print("")
print("Fermi integrals")
yF00 = []
yF05 = []
yF10 = []
yF15 = []
yF20 = []
yG00 = []
yG05 = []
yG10 = []
yG15 = []
yG20 = []
print("{:12}\t{:12}\t{:12}\t{:12}\t{:12}\t{:12}\t{:12}".format("x=EF/kBT", "F0", "F1/2", "F1", "F3/2", "F2" "exp(x)", "Gamma(1)"))
for i in range(nx_basic):
x = xx[i]
for j in range(3):
eta = j / 2.0
y0 = FermiIntegral(x, j)
y1 = FermiIntegral_x2(x, j)
validate_error(y0, y1, 2.0e-6, "Error in basic() for F{}: ".format(eta))
yF00.append(FermiIntegral(x, 0.0))
yF05.append(FermiIntegral(x, 0.5))
yF10.append(FermiIntegral(x, 1.0))
yF15.append(FermiIntegral(x, 1.5))
yF20.append(FermiIntegral(x, 2.0))
yG00.append(exp(x) * gamma(1.0))
yG05.append(exp(x) * gamma(1.5))
yG10.append(exp(x) * gamma(2.0))
yG15.append(exp(x) * gamma(2.5))
yG20.append(exp(x) * gamma(3.0))
if i % 10 == 0:
print("{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}"
.format(x, yF00[i], yF05[i], yF10[i], yF15[i], yF20[i], yG00[i]))
print("")
print("Save data to [{}]".format(outxlsfFDfile))
tkVariousData().to_excel(outxlsfFDfile, ["x=E/kBT", "fFD(exact)", "fFD(approx)", "fFDh(exact)", "-df/dx"],
[xx, yfFD, yfFDapprox, yfFDh, ymdfdx])
print("Save data to [{}]".format(outxlsFjfile))
tkVariousData().to_excel(outxlsFjfile, ["x=EF/kBT", "F0", "F1/2", "F1", "F3/2", "F2", "exp(x)*G(1)",
"exp(x)*G(1.5)", "exp(x)*G(2)", "exp(x)*G(2.5)", "exp(x)*G(3)"],
[xx, yF00, yF05, yF10, yF15, yF20, yG00, yG05, yG10, yG15, yG20])
#=============================
# グラフの表示
#=============================
print("")
fig = plt.figure(figsize = figsize_FD)
ax2 = fig.add_subplot(2, 1, 1)
ax3 = fig.add_subplot(2, 1, 2)
ax2.plot(xx, yfFD, label = '$f_{FD,e}(exact)$', linestyle = '-', color = 'black', linewidth = 0.5)
ax2.plot(xx, yfFDh, label = '$f_{FD,h}(exact)$', linestyle = '-', color = 'blue', linewidth = 0.5)
ax2.plot(xx, yfFDapprox, label = '$f_{FD}(approx)$', linestyle = '-', color = 'green', linewidth = 0.5)
ax2.plot(xx, ymdfdx, label = '$-df/dx$', linestyle = 'dashed', color = 'red', linewidth = 0.5)
ax2.set_xlabel("$x=(E-E_F)/k_BT$ (eV)", fontsize = fontsize)
ax2.set_ylabel("$f_{FD}$, $-df/dx$", fontsize = fontsize)
ax2.set_xlim([-10.0, 10.0])
ax2.set_ylim([-0.1, 1.1])
ax2.legend(fontsize = legend_fontsize)
ax2.tick_params(labelsize = fontsize)
ax3.plot(xx, yF00, label = '$F_0$', linestyle = '-', color = 'black', linewidth = 0.5)
ax3.plot(xx, yF05, label = '$F_{1/2}$', linestyle = '-', color = 'red', linewidth = 0.5)
ax3.plot(xx, yF10, label = '$F_1$', linestyle = '-', color = 'blue', linewidth = 0.5)
ax3.plot(xx, yF15, label = '$F_{3/2}$', linestyle = '-', color = 'purple', linewidth = 0.5)
ax3.plot(xx, yF20, label = '$F_2$', linestyle = '-', color = 'green', linewidth = 0.5)
ax3.plot(xx, yG00, label = r'$e^x$$\Gamma(1)$', linestyle = '', marker = 'o', markerfacecolor = 'black', markersize = 1)
ax3.plot(xx, yG05, label = r'$e^x$$\Gamma(3/2)$', linestyle = '', marker = 'o', markerfacecolor = 'red', markersize = 1)
ax3.plot(xx, yG10, label = r'$e^x$$\Gamma(2)$', linestyle = '', marker = 'o', markerfacecolor = 'blue', markersize = 1)
ax3.plot(xx, yG15, label = r'$e^x$$\Gamma(5/2)$', linestyle = '', marker = 'o', markerfacecolor = 'purple', markersize = 1)
ax3.plot(xx, yG20, label = r'$e^x$$\Gamma(3)$', linestyle = '', marker = 'o', markerfacecolor = 'green', markersize = 1)
ax3.set_xlabel("$x=E_F/k_BT$ (eV)", fontsize = fontsize)
ax3.set_ylabel("$F_r$", fontsize = fontsize)
ax3.set_xlim([-10.0, 10.0])
ax3.set_ylim([1.0e-5, 0.1e4])
ax3.set_yscale('log')
ax3.legend(fontsize = legend_fontsize)
ax3.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
app.terminate("", pause = True)
[ドキュメント]
def Hall():
"""
概要:
ホール係数、ホール因子、ホール移動度などのホール関連特性を計算し、プロットする。
詳細説明:
指定された温度 `T0` における、様々なフェルミ準位 `EF` ( `xmin` から `xmax` の範囲)に対する
キャリア濃度 `N`、電気伝導度 `sigma`、移動度 `mu`、ホール因子 `FH`、ホール係数 `RH`
を計算する。結果はExcelファイル (`outxlsPropfile`) に出力し、matplotlibでグラフ表示する。
引数:
なし。
"""
print("")
print("Carrier:")
print(" q={} e".format(mobility.charge))
print(" meff={}me".format(dos.meeff))
dos.NC = meff2NC_FEA(dos.meeff, T0)
dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
if mobility.charge < 0.0:
print(" NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
print(" DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
else:
print(" NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
print(" DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
print(" scattering factor (non-degenerated) r={}".format(mobility.rfac))
print(" mean free path prefactor: l0={} m".format(mobility.l0))
print("")
print(f"Calculat Hall coefficient, Hall factor, and Hall mobility at {T0} K")
xx = []
yEF = []
yNe = []
ysigma = []
ymu = []
ytau = []
yRH0 = []
yRH = []
yFH = []
yNe_Hall = []
ymu_Hall = []
print("{:>8} {:>8} {:>12} {:>12} {:>12} {:>12} {:>12}"
.format("x=EF/kBT", "EF(eV)", "N(cm^-3)", "sigma(S/cm)", "mu(cm2/Vs)", "FH", "RH(m^3/C"))
for i in range(nx):
x = xmin + i * xstep
EF = x * kB * T0 / e
inf = dos.cal_Hall_properteis(T0, EF, mobility, validate_error_str = 'Error in properteis()')
sigma = inf["sigma"]
Ne = inf["n"]
mu = inf["mu"]
tau = inf["<tau>"]
FH = inf["FH"]
RH = inf["RH"]
RH0 = inf["RH0"]
n_Hall = inf["nHall"]
mu_Hall = inf["muHall"]
xx.append(x)
yEF.append(EF)
yNe.append(Ne)
ysigma.append(sigma)
ymu.append(mu)
ytau.append(FH)
yRH0.append(RH0)
yRH.append(RH)
yFH.append(FH)
yNe_Hall.append(n_Hall)
ymu_Hall.append(mu_Hall)
if i % 10 == 0:
print("{:8.3g} {:8.3g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g}"
.format(x, EF, yNe[i], ysigma[i], ymu[i], yFH[i], yRH[i]))
print("")
print("Save data to [{}]".format(outxlsPropfile))
tkVariousData().to_excel(outxlsPropfile, ["x=EF/kBT", "EF(eV)", "N(cm^-3)", "sigma(S/cm)", "mu(cm2/Vs)", "FH", "RH(m^3/C"],
[xx, yEF, yNe, ysigma, ymu, yFH, yRH])
#=============================
# グラフの表示
#=============================
print("")
fig = plt.figure(figsize = figsize)
axtau = fig.add_subplot(3, 3, 1)
axNe = fig.add_subplot(3, 3, 2)
axsigma = fig.add_subplot(3, 3, 3)
axmu = fig.add_subplot(3, 3, 4)
aRH = fig.add_subplot(3, 3, 5)
aFH = fig.add_subplot(3, 3, 6)
aNeRH = fig.add_subplot(3, 3, 7)
aNeFH = fig.add_subplot(3, 3, 8)
axsigma.set_title("$T_0$={} K $m^*$={}$m_e$ r={} $k_l$$_a$$_t$$_t$={} W/m/K".format(T0, dos.meeff, mobility.rfac, klatt))
axtau.plot(xx, ytau, label = '$\\tau(EF)$', linestyle = '-', color = 'red', linewidth = 0.5)
# axtau.plot(xx, ytauavg, label = '<$\\tau$>', linestyle = '-', color = 'blue', linewidth = 0.5)
axtau.set_xlabel("$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
axtau.set_ylabel("$\\tau$ (fs)", fontsize = fontsize)
axtau.set_ylim([0.0, max(ytau) * 1.1])
axtau.legend(fontsize = legend_fontsize, loc = 'best')
axtau.tick_params(labelsize = fontsize)
axNe.plot(xx, yNe, label = '$N_e$', linestyle = '-', color = 'red', linewidth = 1.0)
axNe.plot(xx, yNe_Hall, label = '$N_e$$_{,Hall}$', linestyle = '-', color = 'blue', linewidth = 1.0)
axNe.set_xlabel("$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
axNe.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
axNe.set_yscale('log')
axNe.legend(fontsize = legend_fontsize, loc = 'best')
axNe.tick_params(labelsize = fontsize)
# axsigma.plot(xx, ysigma, label = '$\sigma$', linestyle = '-', color = 'red', linewidth = 1.0)
axsigma.plot(yNe_Hall, ysigma, label = r'$\sigma$', linestyle = '-', color = 'red', linewidth = 1.0)
# axsigma.set_xlabel(r"$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
axsigma.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
axsigma.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
axsigma.set_xscale('log')
axsigma.set_yscale('log')
axsigma.set_xlim([1.0e10, 1.0e23])
axsigma.legend(fontsize = legend_fontsize, loc = 'best')
axsigma.tick_params(labelsize = fontsize)
# axmu.plot(xx, ymu, label = r'$\mu$', linestyle = '-', color = 'red', linewidth = 1.0)
axmu.plot(yNe_Hall, ymu, label = r'$\mu$', linestyle = '-', color = 'red', linewidth = 1.0)
# axmu.plot(xx, ymu_Hall, label = r'$\mu$$_{,Hall}$', linestyle = '-', color = 'blue', linewidth = 1.0)
axmu.plot(yNe_Hall, ymu_Hall, label = r'$\mu$$_{,Hall}$', linestyle = '-', color = 'blue', linewidth = 1.0)
# axmu.set_xlabel(r"$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
axmu.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
axmu.set_ylabel(r"$\mu$ (cm$^2$/V/s)", fontsize = fontsize)
axmu.set_xscale('log')
axmu.set_xlim([1.0e10, 1.0e23])
axmu.legend(fontsize = legend_fontsize, loc = 'best')
axmu.tick_params(labelsize = fontsize)
# aRH.plot(xx, yRH0, label = r'$R_H$$_0$=$1/e/N_e$', linestyle = '-', color = 'red', linewidth = 1.0)
aRH.plot(yNe_Hall, yRH0, label = r'$R_H$$_0$=$1/e/N_e$', linestyle = '-', color = 'red', linewidth = 1.0)
# aRH.plot(xx, yRH, label = r'$R_H$$_{,Hall}$', linestyle = '-', color = 'blue', linewidth = 1.0)
aRH.plot(yNe_Hall, yRH, label = r'$R_H$$_{,Hall}$', linestyle = '-', color = 'blue', linewidth = 1.0)
# aRH.set_xlabel(r"$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
aRH.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
aRH.set_ylabel(r"$R_H$ (cm$^3$/C)", fontsize = fontsize)
aRH.set_xscale('log')
aRH.set_yscale('log')
aRH.set_xlim([1.0e10, 1.0e23])
aRH.legend(fontsize = legend_fontsize, loc = 'best')
aRH.tick_params(labelsize = fontsize)
# aFH.plot(xx, yFH, label = r'$F_{Hall}$', linestyle = '-', color = 'red', linewidth = 1.0)
aFH.plot(yNe_Hall, yFH, label = r'$F_{Hall}$', linestyle = '-', color = 'red', linewidth = 1.0)
# aFH.set_xlabel(r"$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
aFH.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
aFH.set_ylabel("$F_H$", fontsize = fontsize)
aFH.set_xscale('log')
aFH.set_xlim([1.0e10, 1.0e23])
aFH.legend(fontsize = legend_fontsize, loc = 'best')
aFH.tick_params(labelsize = fontsize)
aNeRH.plot(yNe_Hall, yRH, label = r'$R_H$$_{,Hall}$', linestyle = '-', color = 'red', linewidth = 1.0)
aNeRH.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
aNeRH.set_ylabel(r"$R_H$ (cm$^3$/C)", fontsize = fontsize)
aNeRH.set_xscale('log')
aNeRH.set_yscale('log')
aNeRH.set_xlim([1.0e10, 1.0e23])
# aNeRH.set_ylim([min(yRH) * 0.5, min(yRH) * 1.0e3])
# aNeRH.legend(fontsize = legend_fontsize, loc = 'best')
aNeRH.tick_params(labelsize = fontsize)
aNeFH.plot(yNe_Hall, yFH, label = r'$N_e$', linestyle = '-', color = 'red', linewidth = 1.0)
aNeFH.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
aNeFH.set_ylabel(r"$F_H$ (cm$^{-3}$)", fontsize = fontsize)
aNeFH.set_xscale('log')
aNeFH.set_xlim([1.0e10, 1.0e23])
# aNeFH.legend(fontsize = legend_fontsize, loc = 'best')
aNeFH.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
app.terminate("", pause = True)
[ドキュメント]
def properties():
"""
概要:
ゼーベック係数、電気伝導度、移動度、ローレンツ数、ZT値などの熱電特性を計算し、プロットする。
詳細説明:
指定された温度 `T0` における、様々なフェルミ準位 `EF` ( `xmin` から `xmax` の範囲)に対する
熱電特性を計算する。キャリア濃度、導電率、移動度、平均緩和時間、ゼーベック係数
（正確な計算値、非縮退近似、縮退近似）、電子熱伝導率、全熱伝導率、ローレンツ数、
出力因子、ZT値を算出し、Excelファイル (`outxlsPropfile`) に出力し、matplotlibでグラフ表示する。
引数:
なし。
"""
global mobility, dos
print("")
print("mode:", mode)
print("T={} K".format(T0))
print("meff:", dos.meeff)
dos.NC = meff2NC_FEA(dos.meeff, T0)
dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
if mobility.charge < 0.0:
print(" NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
print(" DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
else:
print(" NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
print(" DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
print(" scattering factor in tau(E) prop to E^(r-0.5) (non-degenerated) r={}".format(mobility.rfac))
print("Lorentz number")
print(" Free electron mode: {:12.8g} Wohm/K^2".format(L_FEA))
print("Mean free path prefactor: l0={} m".format(mobility.l0))
print("Lattice thermal conductivity: klatt={} W/m/K".format(klatt))
print("")
print("Properties")
xx = []
xx_tau = []
yEF = []
yNe = []
ysigma = []
ymu = []
ytau = []
ytauavg = []
yS = []
ySndeg = []
ySdeg = []
ykappa = []
ykappatot = []
yL = []
yLapprox = []
yPF = []
yZT = []
print("{:>8} {:>8} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12}"
.format("x=EF/kBT", "EF(eV)", "N(cm^-3)", "sigma(S/cm)", "mu(cm2/Vs)", "<tau>(fs)",
"S(uV/K)", "Snon-deg", "Sdeg", "kappa,e(W/m/K)", "kappa,tot", "L(Wohm/K^2)", "Lapprox", "PF(W/m/K^2)", "ZT"))
for i in range(nx):
x = xmin + i * xstep
EF = x * kB * T0 / e
# print("T=", T0)
sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
= dos.cal_transport_properteis(T0, EF, mobility, klatt, validate_error_str = 'Error in properteis()')
# Relaxation time tau(EF)
if EF <= 0.0:
tau = None
else:
tau = mobility.tau(T0, EF) * 1.0e15 # fs, E in eV
if EF > 0.0:
xx_tau.append(x)
ytau.append(tau)
Sndeg = dos.cal_S_nondegenerated_from_Ne(n, mobility.rfac, charge = mobility.charge) * 1.0e6 # uV/K
Sdeg = dos.cal_S_degenerated_from_Ne(T0, n, mobility.rfac, charge = mobility.charge) * 1.0e6 # uV/K
# Approximated Lorentz number
Lapprox = dos.cal_LorentzNumber_from_S_approximate(S)
xx.append(x)
yEF.append(EF)
yNe.append(n)
ysigma.append(sigma)
ymu.append(mu)
ytauavg.append(tau_avg)
yS.append(S)
ySndeg.append(Sndeg)
ySdeg.append(Sdeg)
ykappa.append(kappa)
ykappatot.append(kappa + klatt)
yL.append(L)
yLapprox.append(Lapprox)
yPF.append(PF)
yZT.append(ZT)
if i % 10 == 0:
print("{:8.3g} {:8.3g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g}"
.format(x, EF, n, sigma, mu, tau_avg, S, Sndeg, Sdeg, kappa, kappa_tot, L, Lapprox, PF, ZT))
print("")
print("Save data to [{}]".format(outxlsPropfile))
tkVariousData().to_excel(outxlsPropfile, ["x=EF/kBT", "EF(eV)", "N(cm^-3)", "sigma(S/cm)", "mu(cm2/Vs)", "<tau>(fs)",
"S(uV/K)", "Snon-deg", "Sdeg", "kappa,e(W/m/K)", "kappa,tot(W/m/K)", "L(Wohm/K^2)", "Lapprox", "PF(uW/cm/K^2)", "ZT"],
[xx, yEF, yNe, ysigma, ymu, ytauavg, yS, ySndeg, ySdeg, ykappa, ykappatot, yL, yLapprox, yPF, yZT])
#=============================
# グラフの表示
#=============================
print("")
maxS = max(yS)
minS = min(yS)
if maxS > 0.0:
maxS *= 1.2
else:
maxS = 0.0
if minS > 0.0:
minS = 0.0
else:
minS *= 2.0
fig = plt.figure(figsize = figsize)
axtau = fig.add_subplot(3, 4, 1)
axNe = fig.add_subplot(3, 4, 2)
axsigma = fig.add_subplot(3, 4, 3)
axmu = fig.add_subplot(3, 4, 4)
axS = fig.add_subplot(3, 4, 5)
axke = fig.add_subplot(3, 4, 6)
axPF = fig.add_subplot(3, 4, 7)
axZT = fig.add_subplot(3, 4, 8)
axL1 = fig.add_subplot(3, 4, 9)
axL2 = fig.add_subplot(3, 4, 10)
axLS = fig.add_subplot(3, 4, 11)
axsigma.set_title(r"$T_0$={} K $m^*$={}$m_e$ r={} $k_l$$_a$$_t$$_t$={} W/m/K".format(T0, dos.meeff, mobility.rfac, klatt))
axtau.plot(xx_tau, ytau, label = '$\\tau(EF)$', linestyle = '-', color = 'red', linewidth = 0.5)
axtau.plot(xx, ytauavg, label = '<$\\tau$>', linestyle = '-', color = 'blue', linewidth = 0.5)
axtau.set_xlabel("$E$ (eV)", fontsize = fontsize)
axtau.set_ylabel("$\\tau$ (fs)", fontsize = fontsize)
axtau.set_ylim([0.0, max(ytau) * 1.1])
axtau.legend(fontsize = legend_fontsize, loc = 'best')
axtau.tick_params(labelsize = fontsize)
axNe.plot(xx, yNe, label = '$N_e$', linestyle = '-', color = 'red', linewidth = 1.0)
axNe.set_xlabel("$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
axNe.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
axNe.set_yscale('log')
axNe.legend(fontsize = legend_fontsize, loc = 'best')
axNe.tick_params(labelsize = fontsize)
# axsigma.plot(xx, ysigma, label = r'$\sigma$', linestyle = '-', color = 'red', linewidth = 1.0)
axsigma.plot(yNe, ysigma, label = r'$\sigma$', linestyle = '-', color = 'red', linewidth = 1.0)
# axsigma.set_xlabel(r"$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
axsigma.set_xlabel(r"$N_e$ (cm$^{-3}$)", fontsize = fontsize)
axsigma.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
axsigma.set_xscale('log')
axsigma.set_yscale('log')
axsigma.legend(fontsize = legend_fontsize, loc = 'best')
axsigma.tick_params(labelsize = fontsize)
# axmu.plot(xx, ymu, label = '$\\mu$', linestyle = '-', color = 'red', linewidth = 1.0)
axmu.plot(yNe, ymu, label = '$\\mu$', linestyle = '-', color = 'red', linewidth = 1.0)
# axmu.set_xlabel("$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
axmu.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
axmu.set_ylabel("$\\mu$ (cm$^2$/V/s)", fontsize = fontsize)
axmu.set_xscale('log')
axmu.legend(fontsize = legend_fontsize, loc = 'best')
axmu.tick_params(labelsize = fontsize)
axke.plot(yNe, ykappa, label = '$\\kappa$$_e$', linestyle = '-', color = 'red', linewidth = 1.0)
axke.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
# axke.set_xlabel("$x=(E_F-E_C)/k_BT$", fontsize = fontsize)
axke.set_ylabel("$\\kappa$$_e$ (W/m/K)", fontsize = fontsize)
axke.set_xscale('log')
axke.set_yscale('log')
axke.legend(fontsize = legend_fontsize, loc = 'best')
axke.tick_params(labelsize = fontsize)
axS.plot(yNe, yS, label = '$S$', linestyle = '-', color = 'black', linewidth = 1.0)
axS.plot(yNe, ySndeg, label = '$S_{non-deg}$', linestyle = 'dashed', color = 'red', linewidth = 0.5)
axS.plot(yNe, ySdeg, label = '$S_{deg}$', linestyle = 'dashed', color = 'blue', linewidth = 0.5)
axS.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
axS.set_ylabel(r"$S$ ($\mu$V/K)", fontsize = fontsize)
axS.set_xscale('log')
axS.set_ylim([minS, maxS])
axS.legend(fontsize = legend_fontsize, loc = 'best')
axS.tick_params(labelsize = fontsize)
axPF.plot(yNe, yPF, label = '$PF$', linestyle = '-', color = 'black', linewidth = 1.0)
axPF.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
axPF.set_ylabel(r"$PF$ ($\mu$W/m/K$^2$)", fontsize = fontsize)
axPF.set_xscale('log')
axPF.legend(fontsize = legend_fontsize, loc = 'best')
axPF.tick_params(labelsize = fontsize)
axZT.plot(yNe, yZT, label = '$ZT$', linestyle = '-', color = 'black', linewidth = 1.0)
axZT.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
axZT.set_ylabel("$ZT$", fontsize = fontsize)
axZT.set_xscale('log')
axZT.legend(fontsize = legend_fontsize, loc = 'best')
axZT.tick_params(labelsize = fontsize)
# axL1.plot(xx, yL, label = 'L(cal)', linestyle = '-', linewidth = 0.5)
axL1.plot(yNe, yL, label = 'L(cal)', linestyle = '-', linewidth = 0.5)
axL1.plot([min(xx), max(xx)], [L_FEA, L_FEA], label = 'L(FEA)', linestyle = '-', color = 'red', linewidth = 0.5)
# axL1.set_xlabel("$x=(E_F-E_C)/k_B/T$", fontsize = fontsize)
axL1.set_xlabel(r"$N_e$ (cm$^{-3}$)", fontsize = fontsize)
axL1.set_ylabel(r"$L=\kappa_e/\sigma/T$ (W$\Omega$/K$^2$)", fontsize = fontsize)
axL1.set_ylim([0.0, 3.2e-8])
axL1.legend(fontsize = legend_fontsize, loc = 'best')
axL1.tick_params(labelsize = fontsize)
axL2.plot(yNe, yL, label = 'L(cal)', linestyle = '-', linewidth = 0.5)
axL2.plot([min(yNe), max(yNe)], [L_FEA, L_FEA], label = 'L(FEA)', linestyle = '-', color = 'red', linewidth = 0.5)
axL2.set_xlabel(r"$N_e$ (cm$^{-3}$)", fontsize = fontsize)
axL2.set_ylabel(r"$L=\kappa_e/\sigma/T$ (W$\Omega$/K$^2$)", fontsize = fontsize)
axL2.set_xscale('log')
axL2.set_xlim([1.0e10, 1.0e23])
axL2.set_ylim([0.0, 3.2e-8])
axL2.legend(fontsize = legend_fontsize, loc = 'best')
axL2.tick_params(labelsize = fontsize)
axLS.plot(yS, yL, label = 'L(cal)', linestyle = '-', color = 'red', linewidth = 1.0)
axLS.plot(yS, yLapprox, label = 'L(approx, r=0)', linestyle = '-', color = 'blue', linewidth = 0.5)
axLS.plot([min(yS), max(yS)], [L_FEA, L_FEA], label = 'L(FEA)', linestyle = '-', color = 'red', linewidth = 0.5)
axLS.set_xlabel(r"$S$ ($\mu$V/K)", fontsize = fontsize)
axLS.set_ylabel(r"$L=\kappa_e/\sigma/T$ (W$\Omega$/K$^2$)", fontsize = fontsize)
axLS.set_ylim([0.0, 3.2e-8])
axLS.legend(fontsize = legend_fontsize, loc = 'best')
axLS.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
app.terminate("", pause = True)
[ドキュメント]
def recover_parameters(ais, optid, aidef0):
"""
概要:
最適化されたサブセットのパラメータから、元のパラメータリストを再構築する。
詳細説明:
`optid` リストに基づいて、最適化対象ではないパラメータをデフォルト値 (`aidef0`) から取得し、
最適化されたパラメータ (`ais`) と組み合わせて完全なパラメータリストを作成する。
これにより、一部のパラメータのみを最適化した場合でも、全てのパラメータを正しく扱うことができる。
引数:
:param ais: (list) 最適化されたパラメータのサブセット。
:param optid: (list) 各パラメータが最適化対象かどうかを示すIDのリスト (1: 対象, 0: 非対象)。
:param aidef0: (list) デフォルトの全パラメータリスト。
戻り値:
list: 再構築された全パラメータのリスト。
"""
aidef = list(aidef0).copy()
ai = []
c = 0
for i in range(len(optid)):
if optid[i] == 1:
ai.append(ais[c])
c += 1
else:
ai.append(aidef[i])
return ai
[ドキュメント]
def choose_parameters(ais0, optid):
"""
概要:
全パラメータリストから最適化対象のパラメータのみを選択する。
詳細説明:
`optid` リストに基づいて、最適化対象 (`optid` が1) のパラメータのみを抽出して
新しいリストを生成する。これは、`scipy.optimize.minimize` に渡すパラメータを準備するために使用される。
引数:
:param ais0: (list) 全パラメータのリスト。
:param optid: (list) 各パラメータが最適化対象かどうかを示すIDのリスト (1: 対象, 0: 非対象)。
戻り値:
list: 最適化対象のパラメータのみを含むリスト。
"""
ais = list(ais0).copy()
ai = []
c = 0
for i in range(len(optid)):
if optid[i] == 1:
ai.append(ais[i])
return ai
[ドキュメント]
def calS2_sigma(ai):
"""
概要:
実測の電気伝導度と計算値の二乗誤差（`S2`）を計算する。
詳細説明:
与えられた `mobility.l0` の値（`ai`の要素）を用いて、各データポイントにおける
計算された電気伝導度と実測値の差の二乗和を計算する。
`Nmin_fit`, `Nmax_fit`, `sigmamin_fit`, `sigmamax_fit` の範囲外のデータはスキップされる。
計算された二乗誤差はデータ点数で正規化される。
引数:
:param ai: (list) 最適化するパラメータのリスト。現在、`mobility.l0`のみが対象。
戻り値:
float: 電気伝導度の二乗誤差の平均値。
"""
global mu, dos, ysigma
aiorg = parameter_list()
ainew = [ai[0]]
if ainew[0] <= 0.0:
ainew[0] = 1.0e-4
set_parameters([mobility.meff, ainew[0]])
S2 = 0.0
ndata = len(ysigma)
eps = 1.0e-300
for i in range(ndata):
S = yS[i]
sigma = ysigma[i]
n = yN[i]
mu = ymu[i]
if Nmin_fit is not None and n < Nmin_fit:
continue
if Nmax_fit is not None and Nmax_fit < n:
continue
if sigmamin_fit is not None and sigma < sigmamin_fit:
continue
if sigmamax_fit is not None and sigmamax_fit < sigma:
continue
EFfin, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100)
if not flag:
terminate("Error in callback(): EF calculation did not converge. diffEF={} eps={} in {} iteration"
.format(diffEF, epsEF, maxiter))
sigmafin, nfin, mufin, tau_avg, Sfin, kappa, kappa_tot, L, PF, ZT, inf \
= dos.cal_transport_properteis(T0, EFfin, mobility, klatt, validate_error_str = 'Error in fit()')
# print("m,l0=", mobility.meff, mobility.l0)
# print(" sigma=", i, ysigma[i], sigmafin)
dlogsigma = ysigma[i] - sigmafin
# dlogsigma = log(ysigma[i]+eps) - log(sigmafin+eps)
S2 += dlogsigma * dlogsigma
S2 /= (ndata - 1)
set_parameters(aiorg)
# print("S2=", S2)
return S2
[ドキュメント]
def calS2_S(ai):
"""
概要:
実測のゼーベック係数と計算値の二乗誤差（`S2`）を計算する。
詳細説明:
与えられた `dos.meeff` の値（`ai`の要素）を用いて、各データポイントにおける
計算されたゼーベック係数と実測値の差の二乗和を計算する。
`Nmin_fit`, `Nmax_fit`, `sigmamin_fit`, `sigmamax_fit` の範囲外のデータはスキップされる。
計算された二乗誤差はデータ点数で正規化される。
引数:
:param ai: (list) 最適化するパラメータのリスト。現在、`dos.meeff`のみが対象。
戻り値:
float: ゼーベック係数の二乗誤差の平均値。
"""
global mu, dos, ysigma
global method, tol, maxiter, h_diff
global Nmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fit
aiorg = parameter_list()
ainew = [ai[0]]
if ainew[0] <= 0.0:
ainew[0] = 1.0e-4
set_parameters(ainew)
# print_parameters()
# print("m=", dos.meeff, dos.NC, dos.DC0)
S2 = 0.0
ndata = len(ysigma)
eps = 1.0e-300
for i in range(ndata):
S = yS[i]
sigma = ysigma[i]
n = yN[i]
mu = ymu[i]
if Nmin_fit is not None and n < Nmin_fit:
continue
if Nmax_fit is not None and Nmax_fit < n:
continue
if sigmamin_fit is not None and sigma < sigmamin_fit:
continue
if sigmamax_fit is not None and sigmamax_fit < sigma:
continue
EFfin, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100)
if not flag:
app.terminate("Error in callback(): EF calculation did not converge. diffEF={} eps={} in {} iteration"
.format(diffEF, epsEF, maxiter), usage = usage, pause = True)
sigmafin, nfin, mufin, tau_avg, Sfin, kappa, kappa_tot, L, PF, ZT, inf \
= dos.cal_transport_properteis(T0, EFfin, mobility, klatt, validate_error_str = 'Error in fit()')
dlogS = yS[i] - Sfin
"""
if yS[i] <= 0.0:
dlogS = log(-yS[i]+eps) - log(-Sfin+eps)
else:
dlogS = log(yS[i]+eps) - log(Sfin+eps)
"""
S2 += dlogS * dlogS
S2 /= (ndata - 1)
set_parameters(aiorg)
return S2
# First derivatives to be used e.g. for cg method
# Approximate by forward difference method with the delta h = h_diff
[ドキュメント]
def diff1_sigma(ai):
"""
概要:
`calS2_sigma` 関数の数値微分（1階導関数）を計算する。
詳細説明:
前方差分法を用いて、`calS2_sigma` 関数に対する各パラメータの勾配を計算する。
微小な変化量 `h_diff` を使用し、各パラメータを個別に微小量だけ変化させて
目的関数の変化を評価する。
引数:
:param ai: (list) パラメータのリスト。
戻り値:
:returns: (numpy.ndarray) 各パラメータに対する1階導関数の配列。
"""
n = len(ai)
f0 = calS2_sigma(ai)
df = np.empty(n)
for i in range(n):
aii = ai
aii[i] = ai[i] + h_diff
df[i] = (calS2_sigma(aii) - f0) / h_diff
return df
[ドキュメント]
def diff1_S(ai):
"""
概要:
`calS2_S` 関数の数値微分（1階導関数）を計算する。
詳細説明:
前方差分法を用いて、`calS2_S` 関数に対する各パラメータの勾配を計算する。
微小な変化量 `h_diff` を使用し、各パラメータを個別に微小量だけ変化させて
目的関数の変化を評価する。
引数:
:param ai: (list) パラメータのリスト。
戻り値:
:returns: (numpy.ndarray) 各パラメータに対する1階導関数の配列。
"""
n = len(ai)
f0 = calS2_S(ai)
df = np.empty(n)
for i in range(n):
aii = ai
aii[i] = ai[i] + h_diff
df[i] = (calS2_S(aii) - f0) / h_diff
return df
[ドキュメント]
def plot(fig, yn, ysigma, ymu, yS, ysigmaini, ymuini, ySini,
ysigmafin = None, ymufin = None, ySfin = None):
"""
概要:
フィット結果をグラフ表示する。
詳細説明:
ピサレンコプロット (キャリア濃度 vs ゼーベック係数) と
ヨンカープロット (電気伝導度 vs ゼーベック係数) を描画する。
実測値 (`yS`, `ysigma`, `yN`) と初期計算値 (`ySini`, `ysigmaini`)、
最終的なフィット結果 (`ySfin`, `ysigmafin`) を比較表示する。
グラフはリアルタイムで更新されるように設計されている。
引数:
:param fig: (matplotlib.figure.Figure) 描画対象のFigureオブジェクト。
:param yn: (list) キャリア濃度データ。
:param ysigma: (list) 電気伝導度データ。
:param ymu: (list) 移動度データ。
:param yS: (list) ゼーベック係数データ（実測値）。
:param ysigmaini: (list) 初期計算の電気伝導度データ。
:param ymuini: (list) 初期計算の移動度データ。
:param ySini: (list) 初期計算のゼーベック係数データ。
:param ysigmafin: (list, optional) 最終フィットの電気伝導度データ。デフォルトは None。
:param ymufin: (list, optional) 最終フィットの移動度データ。デフォルトは None。
:param ySfin: (list, optional) 最終フィットのゼーベック係数データ。デフォルトは None。
"""
global plt
global graphupdateinterval
plt.clf()
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)
plt.title(infile, fontsize = fontsize)
Srange = [min(yS) * 0.9, max(yS) * 1.1]
yn2, ySi2 = sort_lists([yn, ySini])
ins1 = ax1.plot(yn, yS, label = 'Pisarenko(obs)', linestyle = 'none', marker = 'o', markerfacecolor = 'blue', markeredgecolor = 'blue')
ins2 = ax1.plot(yn2, ySi2, label = 'Pisarenko(ini)', linestyle = '-', color = 'blue', linewidth = 0.5)
# ins2 = ax1.plot(yn, ySini, label = 'Pisarenko(ini)', linestyle = '-', color = 'blue', linewidth = 0.5)
if ysigmafin:
yn2, ySf2 = sort_lists([yn, ySini])
ins3 = ax1.plot(yn2, ySf2, label = 'Pisarenko(fin)', linestyle = '-', color = 'red')
# ins3 = ax1.plot(yn, ySfin, label = 'Pisarenko(fin)', linestyle = '-', color = 'red')
ax1.set_xscale('log')
ax1.set_ylim([0.0, Srange[1]])
ax1.set_xlabel(r'$N_e$ (cm$^{-3}$)', fontsize = fontsize)
ax1.set_ylabel(r'S ($\mu$V/K)', fontsize = fontsize)
ax1.legend(fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
ax1.tick_params(labelsize = fontsize)
ysi, ySi = sort_lists([ysigmaini, ySini])
ins1 = ax2.plot(ysigma, yS, label = 'Jonker(obs)', linestyle = 'none', marker = 'o', markerfacecolor = 'blue', markeredgecolor = 'blue')
ins2 = ax2.plot(ysi, ySi, label = 'Jonker(ini)', linestyle = '-', color = 'blue', linewidth = 0.5)
# ins2 = ax2.plot(ysigmaini, ySini, label = 'Jonker(ini)', linestyle = '-', color = 'blue', linewidth = 0.5)
if ysigmafin:
ysf, ySf = sort_lists([ysigmafin, ySfin])
ins3 = ax2.plot(ysf, ySf, label = 'Jonker(fin)', linestyle = '-', color = 'red')
# ins3 = ax2.plot(ysigmafin, ySfin, label = 'Jonker(fin)', linestyle = '-', color = 'red')
ax2.set_xscale('log')
ax2.set_ylim([0.0, Srange[1]])
ax2.set_xlabel(r'$\sigma$ (S/cm)', fontsize = fontsize)
ax2.set_ylabel(r'S ($\mu$V/K)', fontsize = fontsize)
ax2.legend(fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
ax2.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.01)
# Callback function for scipy.optimize.minimize()
# Print variables every iteration, and update graph for every graphupdateinterval iterationsおｔ
iter = 0
[ドキュメント]
def callback(fig, S2):
"""
概要:
`scipy.optimize.minimize` のコールバック関数。最適化の進行状況をコンソール出力し、グラフを更新する。
詳細説明:
`outputinterval` ごとに現在のパラメータとS2値をコンソールに出力し、
`graphupdateinterval` ごとに現在のフィット結果をプロットで更新する。
各データポイントのEFを再計算し、それに基づいて`sigmafin`, `mufin`, `Sfin` を求め、
`plot` 関数に渡してグラフを更新する。
引数:
:param fig: (matplotlib.figure.Figure) 描画対象のFigureオブジェクト。
:param S2: (float) 現在の二乗誤差。
"""
global iter, graphupdateinterval
global outputinterval
global dos, mu
global ysigma
global method, tol, maxiter, h_diff
global Nmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fit
if iter % outputinterval == 0:
print("callback {:04d}: {:10.4g} {:10.4g} S2={:10.6g}".format(iter, mobility.meff, mobility.l0, S2))
if iter % graphupdateinterval == 0:
ndata = len(ysigma)
ysigmafin = []
ymufin = []
ySfin = []
for i in range(ndata):
S = yS[i]
sigma = ysigma[i]
n = yN[i]
mu = ymu[i]
if Nmin_fit is not None and n < Nmin_fit:
continue
if Nmax_fit is not None and Nmax_fit < n:
continue
if sigmamin_fit is not None and sigma < sigmamin_fit:
continue
if sigmamax_fit is not None and sigmamax_fit < sigma:
continue
EFfin, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100)
if not flag:
app.terminate("Error in callback(): EF calculation did not converge. diffEF={} eps={} in {} iteration"
.format(diffEF, epsEF, maxiter), usage = usage, pause = True)
sigmafin, nfin, mufin, tau_avg, Sfin, kappa, kappa_tot, L, PF, ZT, inf \
= dos.cal_transport_properteis(T0, EFfin, mobility, klatt, validate_error_str = 'Error in fit()')
ysigmafin.append(sigmafin)
ymufin.append(mufin)
ySfin.append(Sfin)
plot(fig, yN, ysigma, ymu, yS, ysigmaini, ymuini, ySini, ysigmafin, ymufin, ySfin)
iter += 1
[ドキュメント]
def callback_S(fig, xk):
"""
概要:
ゼーベック係数の最適化におけるコールバック関数。
詳細説明:
`scipy.optimize.minimize` が`dos.meeff`の最適化を実行する際に呼び出される。
`calS2_S` を呼び出して現在の `xk` (有効質量)に対するS2を計算し、
`set_parameters` で`dos.meeff`を更新後、汎用コールバック関数 `callback` を呼び出して
出力とグラフ更新を行う。
引数:
:param fig: (matplotlib.figure.Figure) 描画対象のFigureオブジェクト。
:param xk: (numpy.ndarray) 現在の最適化パラメータ（`dos.meeff`）。
"""
S2 = calS2_S(xk)
ainew = [xk[0]]
set_parameters(ainew)
callback(fig, S2)
[ドキュメント]
def callback_sigma(fig, xk):
"""
概要:
電気伝導度の最適化におけるコールバック関数。
詳細説明:
`scipy.optimize.minimize` が`mobility.l0`の最適化を実行する際に呼び出される。
`calS2_sigma` を呼び出して現在の `xk` (平均自由行程のプレファクター)に対するS2を計算し、
`set_parameters` で`mobility.l0`を更新後、汎用コールバック関数 `callback` を呼び出して
出力とグラフ更新を行う。
引数:
:param fig: (matplotlib.figure.Figure) 描画対象のFigureオブジェクト。
:param xk: (numpy.ndarray) 現在の最適化パラメータ（`mobility.l0`）。
"""
global iter, graphupdateinterval
global outputinterval
global dos, mu
global ysigma
S2 = calS2_sigma(xk)
ainew = [mobility.meff, xk[0]]
set_parameters(ainew)
callback(fig, S2)
[ドキュメント]
def construct_lists(label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu):
"""
概要:
読み込んだデータを解析し、不足しているデータを補完する。
詳細説明:
`read_datafile` から返されたリストを整理し、もし `label_sample`, `label_sigma`,
`label_N`, `label_mu` が None の場合、デフォルト値や他のデータからの計算でリストを
生成・補完する。例えば、`yN` が None で `ysigma` と `ymu` が存在する場合、
`yN` は `ysigma / (e * ymu)` から計算される。
引数:
:param label_sample: (str) サンプル名ラベル。
:param xsample: (list) サンプルデータ。
:param label_S: (str) ゼーベック係数ラベル。
:param yS: (list) ゼーベック係数データ。
:param label_sigma: (str) 電気伝導度ラベル。
:param ysigma: (list) 電気伝導度データ。
:param label_N: (str) キャリア濃度ラベル。
:param yN: (list) キャリア濃度データ。
:param label_mu: (str) 移動度ラベル。
:param ymu: (list) 移動度データ。
戻り値:
tuple: 補完されたラベルとデータリストのタプル。
(label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu)
"""
ndata = len(yS)
if label_sample is None:
label_sample = 'sample'
xsample = []
for i in range(ndata):
xsample.append(i + 1)
if label_sigma is None:
ysigma = []
for i in range(ndata):
ysigma.append(None)
if label_N is None:
yN = []
for i in range(ndata):
if ymu is not None:
yN.append(None)
else:
yN.append(ysigma[i] / e / ymu[i])
if label_mu is None:
ymu = []
for i in range(ndata):
if yN[i] is not None:
ymu.append(None)
else:
ymu.append(ysigma[i] / e / yN[i])
return label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu
[ドキュメント]
def fit():
"""
概要:
実測のゼーベック係数データと電気伝導度データにフィッティングを行う。
詳細説明:
入力ファイルからデータを読み込み、指定された温度`T0`と散乱因子`rfac`で、
有効質量`meff`と平均自由行程の係数`l0`をフィッティングパラメータとして最適化する。
`scipy.optimize.minimize` を使用し、`calS2_S` (ゼーベック係数) を目的関数として `meff` を、
`calS2_sigma` (電気伝導度) を目的関数として `l0` をそれぞれ最適化する。
最適化の過程はグラフとコンソールに表示される。
最終的なフィット結果とパラメータはExcelファイルとINIファイルに保存される。
引数:
なし。
"""
global fig
global label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu
global ysigmaini, ymuini, ySini
global method, tol, maxiter, h_diff
global Nmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fit
print("infile :", infile)
print(" T_label : ", T_label)
print(" S_label : ", S_label)
print(" n_label : ", n_label)
print(" mu_label : ", mu_label)
print(" sigma_label: ", sigma_label)
print("")
print(f"T0: {T0} K")
print("Carrier:")
print(" q={} e".format(mobility.charge))
print(" meff={}me".format(dos.meeff))
dos.NC = meff2NC_FEA(dos.meeff, T0)
dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
print(" NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
print(" DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
print(" scattering factor (non-degenerated) r={}".format(mobility.rfac))
print(" mean free path prefactor: l0={} m".format(mobility.l0))
print("")
print( "Fitting condition")
print(f" method : {method}")
print(f" tol : {tol}")
print(f" maxiter: {maxiter}")
print(f" h_diff: {h_diff}")
print(f" outputinterval: {outputinterval}")
print(f" N range : {Nmin_fit} - {Nmax_fit} cm^-3")
print(f" sigma range: {sigmamin_fit} - {sigmamax_fit} S/cm")
def set_default_range(var):
if var == '' or var == '*':
return None
return pfloat(var, defval = None)
Nmin_fit = set_default_range(Nmin_fit)
Nmax_fit = set_default_range(Nmax_fit)
sigmamin_fit = set_default_range(sigmamin_fit)
sigmamax_fit = set_default_range(sigmamax_fit)
if '***' in method:
app.terminate(f"***Error: Choose method", pause = True)
if '***' in S_label:
app.terminate(f"***Error: Choose S_label", pause = True)
if '***' in n_label:
app.terminate(f"***Error: Choose n_label", pause = True)
if '***' in mu_label:
app.terminate(f"***Error: Choose mu_label", pause = True)
if 'uV' in S_label or 'micro' in S_label or 'mV' in S_label:
app.terminate(f"***Error: The unit of S must be 'V/K' but given by [{S_label}]", pause = True)
print("")
print("Read S data from {}".format(infile))
datafile = tkVariousData(infile)
labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = usage)
label_S, yS = datafile.FindDataArray(S_label, flag = 'i')
# Convert V/K to uV/K
for i in range(len(yS)):
yS[i] *= 1.0e6
label_sigma, ysigma = datafile.FindDataArray(sigma_label, flag = 'i')
label_N, yN = datafile.FindDataArray(n_label, flag = 'i')
label_mu, ymu = datafile.FindDataArray(mu_label, flag = 'i')
label_sample, xsample = labels[0], datalist[0]
ndata = len(yS)
print("ndata(all): ", ndata)
ndata_fit = 0
print("data to be fitted")
print(f"{'sample':15} {'n(cm-3)':12} {'mu(cm2/Vs)':10} {'sigma(S/cm)':10} {'S(uV/K)':10}")
for i in range(ndata):
sigma = ysigma[i]
n = yN[i]
if Nmin_fit is not None and n < Nmin_fit:
continue
if Nmax_fit is not None and Nmax_fit < n:
continue
if sigmamin_fit is not None and sigma < sigmamin_fit:
continue
if sigmamax_fit is not None and sigmamax_fit < sigma:
continue
print(f"{xsample[i]:15} {n:12.4g} {ymu[i]:10.4g} {sigma:10.4g} {yS[i]:10.4g}")
print("")
print("Calculate initial values")
print_parameters()
yEFini = []
ysigmaini = []
ymuini = []
ySini = []
print("{:>8} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12}"
.format("EF(eV)", "Ne(cm-3)", "sigma(S/cm)", "mu(cm2/Vs)", "S(uV/K)", "Ne,ini", "sigma,ini", "mu,ini", "S,ini"))
for i in range(ndata):
sample = xsample[i]
S = yS[i]
sigma = ysigma[i]
n = yN[i]
mu = ymu[i]
EFini, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100)
if not flag:
app.terminate("Error in fit(): EF calculation did not converge. diffEF={} eps={} in {} iteration".format(diffEF, epsEF, maxiter),
usage = usage, pause = True)
# x = (EF - doc.EC) * e / kB / T0
sigmaini, nini, muini, tau_avg, Sini, kappa, kappa_tot, L, PF, ZT, inf = \
dos.cal_transport_properteis(T0, EFini, mobility, klatt, validate_error_str = 'Error in fit()')
yEFini.append(EFini)
ysigmaini.append(sigmaini)
ymuini.append(muini)
ySini.append(Sini)
print("{:8.3g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g}"
.format(EFini, n, sigma, mu, S, nini, sigmaini, muini, Sini))
#=============================
# グラフの表示
#=============================
print("")
print("plot")
fig = plt.figure(figsize = figsize_sim)
plot(fig, yN, ysigma, ymu, yS, ysigmaini, ymuini, ySini)
#=============================
# Optimization
#=============================
print("")
print("Nonlinear least-squares fitting for effective mass by ", method)
print(" tol=", tol)
ai = [mobility.meff]
ret = minimize(calS2_S, ai, method = method, jac = diff1_S, tol = tol, callback = lambda xk: callback_S(fig, xk),
options = {'maxiter':maxiter, "disp":True})
if method == 'nelder-mead':
simplex = ret['final_simplex']
ai = simplex[0][0]
fmin = ret['fun']
aifin = [ai[0], mobility.l0] #ai0[1]]
set_parameters(aifin)
print("Optimized at S2={:12.6g}:".format(fmin))
print_parameters(parameter_list())
print("")
print("Nonlinear least-squares fitting for l0 by ", method)
print(" tol=", tol)
ai_sigma = [mobility.l0]
print("meff(ini)=", mobility.meff)
# ret = minimize(calS2_S, ai_sigma, method = method, jac = diff1_sigma, tol = tol, callback = lambda xk: callback_sigma(fig, xk),
ret = minimize(calS2_sigma, ai_sigma, method = method, jac = diff1_sigma, tol = tol, callback = lambda xk: callback_sigma(fig, xk),
options = {'maxiter':maxiter, "disp":True})
if method == 'nelder-mead':
simplex = ret['final_simplex']
ai_sigma = simplex[0][0]
fmin = ret['fun']
aifin = [mobility.meff, ai_sigma[0]]
set_parameters(aifin)
print("Optimized at S2={:12.6g}:".format(fmin))
# print_parameters()
print_parameters(parameter_list())
print("")
print("Calculate final values")
print_parameters(parameter_list())
yEFfin = []
ysigmafin = []
ySfin = []
print("{:>8} {:>12} {:>12} {:>12} {:>12} {:>12}"
.format("EF(eV)", "Ne(cm-3)", "sigma(S/cm)", "mu(cm2/Vs)", "S(uV/K)", "S,fin"))
for i in range(ndata):
sample = xsample[i]
S = yS[i]
sigma = ysigma[i]
n = yN[i]
mu = ymu[i]
EFfin, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100)
if not flag:
app.terminate("Error in fit(): EF calculation did not converge. diffEF={} eps={} in {} iteration".format(diffEF, epsEF, maxiter),
usage = usage, pause = True)
sigmafin, nfin, mufin, tau_avg, Sfin, kappa, kappa_tot, L, PF, ZT, inf \
= dos.cal_transport_properteis(T0, EFfin, mobility, klatt, validate_error_str = 'Error in fit()')
mufin = sigmafin / sigma / mufin
sigmafin = e * nfin * mufin
yEFfin.append(EFfin)
ysigmafin.append(sigmafin)
ySfin.append(Sfin)
print("{:8.3g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g}"
.format(EFfin, n, sigma, mu, S, Sfin))
print("")
print("Save final parameters to [{}]".format(parameterfile))
save_parameterfile(ai = aifin, S2 = fmin)
print("Save fitting data to [{}]".format(outxlsfile))
tkVariousData().to_excel(outxlsfile, ["EF(eV)", "Ne(cm-3)", "sigma(S/cm)", "mu(cm2/Vs)", "S(uV/K)", "S,fin"],
[yEFfin, yN, ysigma, ymu, yS, ySfin])
app.terminate("", pause = True)
[ドキュメント]
def sim():
"""
概要:
入力ファイルからデータを読み込み、熱電特性をシミュレーションし、
ヨンカープロットやピサレンコプロットを表示する。
詳細説明:
実測のゼーベック係数、電気伝導度、キャリア濃度、ホール移動度から加重移動度を計算し、
さらに指定された有効質量と散乱因子を用いて理論的な熱電特性曲線を生成する。
実測データと理論曲線を比較してグラフ表示する。
計算結果はExcelファイル (`outxlsfile`) に保存される。
引数:
なし。
"""
global label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu
global Nmin, Nmax
print("")
print("infile :", infile)
print(" T_label : ", T_label)
print(" S_label : ", S_label)
print(" n_label : ", n_label)
print(" mu_label : ", mu_label)
print(" sigma_label: ", sigma_label)
print(f"T0 : {T0} K")
print("Carrier:")
print(" q={} e".format(mobility.charge))
print(" meff={}me".format(dos.meeff))
dos.NC = meff2NC_FEA(dos.meeff, T0)
dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
if mobility.charge < 0.0:
print(" NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
print(" DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
else:
print(" NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
print(" DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
print(" scattering factor (non-degenerated) r={}".format(mobility.rfac))
print(" mean free path prefactor: l0={} m".format(mobility.l0))
if '***' in S_label:
app.terminate(f"***Error: Choose S_label", pause = True)
if '***' in n_label:
app.terminate(f"***Error: Choose n_label", pause = True)
if '***' in mu_label:
app.terminate(f"***Error: Choose mu_label", pause = True)
if 'uV' in S_label or 'micro' in S_label or 'mV' in S_label:
print("")
print(f"***Error: The unit of S must be 'V/K' but given by [{S_label}]")
app.terminate(pause = True)
print("")
print("Read S data from {}".format(infile))
datafile = tkVariousData(infile)
labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = usage)
label_S, yS = datafile.FindDataArray(S_label, flag = 'i')
# Convert V/K to uV/K
for i in range(len(yS)):
yS[i] *= 1.0e6
label_sigma, ysigma = datafile.FindDataArray(sigma_label, flag = 'i')
label_N, yN = datafile.FindDataArray(n_label, flag = 'i')
label_mu, ymu = datafile.FindDataArray(mu_label, flag = 'i')
label_sample, xsample = labels[0], datalist[0]
if T_label == '' or '***' in T_label:
label_T = 'T(K)'
yT = []
for i in range(len(xsample)):
yT.append(T0)
else:
label_T, yT = datafile.FindDataArray(T_label, flag = 'i')
# print("x=", xsample)
ndata = len(yS)
print("ndata(all): ", ndata)
# print("xsample=", label_sample, xsample)
# print("yS=", label_S, yS)
# print("ysigma=", label_sigma, ysigma)
print(f"{'sample':15} {'n(cm-3)':12} {'mu(cm2/Vs)':10} {'sigma(S/cm)':10} {'S(uV/K)':10}")
for i in range(ndata):
sigma = ysigma[i]
n = yN[i]
print(f"{xsample[i]:15} {n:12.4g} {ymu[i]:10.4g} {sigma:10.4g} {yS[i]:10.4g}")
print("")
print("Weighted mobility")
ysigma_w = []
yn_w = []
ymu_w = []
ymu_w0 = []
yEF = []
if yN[0] is not None:
print("{:>10}\t{:>8}\t{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>15}\t{:>15}"
.format("sample", "T(K)", "EF(eV)", "S(uV/K)", "sigma(S/cm)", "mu(cm2/Vs)", "Ne(cm-3)", "sigma_w", "mu_w(cm2/Vs)", "mu_w*(me/m*)^(3/2)", "n_w(cm-3)"))
else:
print("{:>10}\t{:>10}\t{:>10}\t{:>10}".format("sample", "EF(eV)", "S(uV/K)", "sigma(S/cm)"))
for i in range(ndata):
sample = xsample[i]
S = yS[i]
sigma = ysigma[i]
n = yN[i]
mu = ymu[i]
T = yT[i]
print(f"sample={sample} S={S} uV/K sigma={sigma} S/cm")
mu_w, n_w, sign, carriertype = weighted_mobility(sigma / 0.01, S * 1.0e-6, T0)
mu_w *= 1.0 * 1.0e4 # cm2/Vs
mu_w0 = mu_w * pow(mobility.meff, -1.5)
n_w *= 1.0e-6 # cm^-3
sigma_w = e * n_w * mu_w
if n is not None:
EF, diffEF, flag = dos.EF_from_electrondensity(n, T, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100)
else:
EF, diffEF, flag = dos.EF_from_electronSeebeck(S, T, mobility, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100)
if not flag:
app.terminate("Error in sim(): EF calculation did not converge. diffEF={} eps={} in {} iteration".format(diffEF, epsEF, maxiter),
usage = usage, pause = True)
# print("n=", n, EF)
ysigma_w.append(sigma_w)
yn_w.append(n_w)
ymu_w.append(mu_w)
ymu_w0.append(mu_w0)
yEF.append(EF)
if n is not None:
print("{:10.4g}\t{:8.3g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:15.4g}\t{:15.4g}"
.format(sample, T, EF, S, sigma, mu, n, sigma_w, mu_w, mu_w0, n_w))
else:
print("{:10.4g}\t{:8.3g}\t{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, T, EF, S, sigma))
Nmin_data = min(yN)
Nmax_data = max(yN)
if Nmin_data < Nmin:
Nmin = Nmin_data
if Nmax < Nmax_data:
Nmax = Nmax_data
lnNmin = log(Nmin)
lnNmax = log(Nmax)
lnNstep = (lnNmax - lnNmin) / (nN - 1)
print("")
print("N range: {:10.4g} - {:10.4g} cm-3, {:10.4g} step in ln(N), nN={}".format(Nmin, Nmax, lnNstep, nN))
xNsim = []
xsigmasim = []
yScal = []
ySsim_nondeg = []
ySsim_deg = []
print("{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}"
.format("N(cm-3)", "EF(eV)", "S,cal(uV/K)", "S,non-deg(uV/K)", "S,deg(uV/K)", "sigma,cal(S/cm)", "Ne,cal(cm-3)", "mu,cal(cm2/Vs)"))
for iN in range(nN):
lnN = lnNmin + lnNstep * iN
N = exp(lnN)
# print("N=", N, dos.NC)
if N < dos.NC:
EF = kB * T0 * log(N / dos.NC) / e
else:
EF = BMShift_FEA(N, dos.DC0)
sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
= dos.cal_transport_properteis(T0, EF, mobility, klatt, validate_error_str = 'Error in sim()')
Sndeg = dos.cal_S_nondegenerated_from_Ne(n, mobility.rfac, charge = mobility.charge) * 1.0e6 # uV/K
Sdeg = dos.cal_S_degenerated_from_Ne(T0, n, mobility.rfac, charge = mobility.charge) * 1.0e6 # uV/K
xNsim.append(n)
xsigmasim.append(sigma)
yScal.append(S)
ySsim_nondeg.append(Sndeg)
ySsim_deg.append(Sdeg)
print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}"
.format(N, EF, S, Sndeg, Sdeg, sigma, n, mu))
print("")
print("Save parameters to [{}]".format(parameterfile))
save_parameterfile(S2 = '')
print("")
print("Save data to [{}]".format(outxlsfile))
if yN[0] is not None:
tkVariousData().to_excel(outxlsfile, ["sample", 'S(uV/K)', 'sigma(S/cm)', 'mu(cm2/Vs)', 'N(cm-3)', "sigma_w", "mu_w", "mu_w*(me/m*)^(3/2)", "n_w"],
[xsample, yS, ysigma, ymu, yN, ysigma_w, ymu_w, ymu_w0, yn_w])
else:
tkVariousData().to_excel(outxlsfile, ["sample", 'S(uV/K)', 'sigma(S/cm)'], [xsample, yS, ysigma])
#=============================
# グラフの表示
#=============================
print("")
fig = plt.figure(figsize = figsize_sim)
ax1 = fig.add_subplot(2, 2, 1)
ax1b = ax1.twinx()
ax2 = fig.add_subplot(2, 2, 2)
ax4 = fig.add_subplot(2, 2, 3)
ax4b = ax4.twinx()
ax3 = fig.add_subplot(2, 2, 4)
maxS = max(yS)
minS = min(yS)
if maxS > 0.0:
maxS *= 2.0
else:
maxS = 0.0
if minS > 0.0:
minS = 0.0
else:
minS *= 2.0
ins1 = ax1.plot( xsample, yS, label = 'S', linestyle = '-', linewidth = 1.0, color = 'red', marker = 'o', markersize = 5.0)
ins2 = ax1b.plot(xsample, ysigma, label = r'$\sigma$', linestyle = '-', linewidth = 1.0, color = 'blue', marker = '^', markersize = 10.0)
ins3 = ax1b.plot(xsample, ysigma_w, label = r'$\sigma_w$', linestyle = 'dashed', linewidth = 0.5, color = 'green', marker = '+', markersize = 15.0, markerfacecolor = 'green')
ax1.set_xlabel("sample", fontsize = fontsize)
ax1.set_ylabel(r"S (K/$\mu$V)", fontsize = fontsize)
ax1b.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
# ax1.set_ylim([minS, maxS])
# h1a, l1a = ax1.get_legend_handles_labels()
# h1b, l1b = ax1b.get_legend_handles_labels()
# ax1.legend(h1a + h1b, l1a + h1b, fontsize = legend_fontsize, loc = 0)
ins = ins1 + ins2 + ins3
ax1.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
# ax1.legend(loc = 0)
# ax1b.legend(loc = 0)
ax1.tick_params(labelsize = fontsize)
ax1b.tick_params(labelsize = fontsize)
if yN[0] is not None:
ins1 = ax4.plot( xsample, ymu, label = r'$\mu$', linestyle = '-', color = 'red', linewidth = 0.5, marker = 'o', markersize = 10.0)
ins2 = ax4.plot( xsample, ymu_w, label = r'$\mu_w$', linestyle = 'dashed', color = 'red', linewidth = 0.5, marker = '^', markersize = 5.0)
ins3 = ax4.plot( xsample, ymu_w0, label = r'$\mu_w^0$', linestyle = 'dashed', color = 'red', linewidth = 0.5, marker = '+', markersize = 5.0)
ins4 = ax4b.plot(xsample, yN, label = '$N_e$', linestyle = '-', color = 'blue', linewidth = 0.5, marker = 'o', markersize = 10.0)
ins5 = ax4b.plot(xsample, yn_w, label = '$N_{e,w}$', linestyle = 'dashed', color = 'blue', linewidth = 0.5, marker = '^', markersize = 5.0)
ax4.set_xlabel("sample", fontsize = fontsize)
ax4.set_ylabel(r"$\mu$ (cm$^2$/V/s)", fontsize = fontsize)
ax4b.set_ylabel(r"$N_e$ (cm$^{-3}$)", fontsize = fontsize)
# ax4.set_ylim([minS, maxS])
ax4b.set_yscale('log')
ins = ins1 + ins2 + ins3 + ins4 + ins5
ax4.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
ax4.tick_params(labelsize = fontsize)
ax4b.tick_params(labelsize = fontsize)
ax2.plot(ysigma, yS, label = 'Jonker plot', linestyle = '', color = 'red', marker = 'o', markersize = 5.0)
ax2.plot(xsigmasim, yScal, label = 'S,cal', linestyle = '-', color = 'red', linewidth = 1.0)
ax2.plot(xsigmasim, ySsim_nondeg, label = 'S,sim,nondeg', linestyle = 'dashed', color = 'blue', linewidth = 0.5)
ax2.plot(xsigmasim, ySsim_deg, label = 'S,sim,deg', linestyle = 'dashed', color = 'green', linewidth = 0.5)
ax2.set_xlabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
ax2.set_ylabel(r"S ($\mu$V/K)", fontsize = fontsize)
ax2.set_xscale('log')
ax2.set_ylim([minS, maxS])
ax2.legend(fontsize = legend_fontsize, loc = 'best')
ax2.tick_params(labelsize = fontsize)
if yN[0] is not None:
ax3.plot(yN, yS, label = 'Pisarenko plot', linestyle = '', marker = 'o', markersize = 5.0)
ax3.plot(xNsim, yScal, label = 'S,cal', linestyle = '-', color = 'red', linewidth = 1.0)
ax3.plot(xNsim, ySsim_nondeg, label = 'S,sim,nondeg', linestyle = 'dashed', color = 'blue', linewidth = 0.5)
ax3.plot(xNsim, ySsim_deg, label = 'S,sim,deg', linestyle = 'dashed', color = 'green', linewidth = 0.5)
ax3.set_xlabel(r"$N$ (cm$^{-3}$)", fontsize = fontsize)
ax3.set_ylabel(r"S ($\mu$V/K)", fontsize = fontsize)
ax3.set_xscale('log')
ax3.set_ylim([minS, maxS])
ax3.legend(fontsize = legend_fontsize, loc = 'best')
ax3.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
app.terminate("", pause = True)
[ドキュメント]
def init():
"""
概要:
最適化のための初期パラメータファイル（.inファイル）を生成する。
詳細説明:
入力ファイルからデータを読み込み、キャリアのタイプ（電子または正孔）を
ゼーベック係数の符号に基づいて自動的に判断する。
その後、現在のグローバルパラメータをデフォルト値としてINIファイルに保存する。
`S`の単位が'V/K'であることを確認し、必要に応じてエラーを発生させる。
引数:
なし。
"""
print("")
print("Carrier:")
print(" q={} e".format(mobility.charge))
print(" meff={}me".format(dos.meeff))
dos.NC = meff2NC_FEA(dos.meeff, T0)
dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
print(" NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
print(" DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
print(" scattering factor (non-degenerated) r={}".format(mobility.rfac))
print(" mean free path prefactor: l0={} m".format(mobility.l0))
if 'uV' in S_label or 'micro' in S_label or 'mV' in S_label:
print("")
print(f"***Error: The unit of S must be 'V/K' but given by [{S_label}]")
app.terminate(pause = True)
print("")
print("Read S data from {}".format(infile))
label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu \
= read_datafile(infile, usage = usage)
# Convert V/K to uV/K
for i in range(len(yS)):
yS[i] *= 1.0e6
ndata = len(xsample)
print("ndata: ", ndata)
if yN[0] is not None:
print("{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}"
.format("sample", "S(uV/K)", "sigma(S/cm)", "mu(cm2/Vs)", "Ne(cm-3)"))
else:
print("{:>10}\t{:>10}\t{:>10}".format("sample", "S(uV/K)", "sigma(S/cm)"))
for i in range(ndata):
sample = xsample[i]
S = yS[i]
sigma = ysigma[i]
n = yN[i]
mu = ymu[i]
if n is not None:
print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, S, sigma, mu, n))
else:
print("{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, S, sigma))
if yS[0] > 0.0:
mobility.charge = 1.0
elif yS[0] < 0.0:
mobility.charge = -1.0
for i in range(1, ndata):
if yS[i] * mobility.charge < 0.0:
app.terminate("Error in init(): S data include both positive and negative values: S[{}}={}, S[{}]={}".format(0, yS[0], i, yS[i]),
usage = usage, pause = True)
print("")
print("Save parameters to [{}]".format(parameterfile))
save_parameterfile(S2 = '')
[ドキュメント]
def calL():
"""
概要:
実測の電気伝導度、キャリア濃度、温度からローレンツ数と電子熱伝導率を計算する。
詳細説明:
入力ファイルから温度、キャリア濃度、電気伝導度のデータを読み込み、
各データポイントにおけるフェルミ準位を計算する。
そのフェルミ準位を用いて理論的なローレンツ数と電子熱伝導率を導出し、
結果をExcelファイル (`outxlsLfile`) に出力し、グラフ表示する。
引数:
なし。
"""
print("")
print("Carrier:")
print(" q={} e".format(mobility.charge))
print(" meff={}me".format(dos.meeff))
dos.NC = meff2NC_FEA(dos.meeff, T0)
dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
if mobility.charge < 0.0:
print(" NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
print(" DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
else:
print(" NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
print(" DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
print(" scattering factor (non-degenerated) r={}".format(mobility.rfac))
print(" mean free path prefactor: l0={} m".format(mobility.l0))
print("")
print("Read T, Ne, and sigma data from {}".format(infile))
datafile = tkVariousData(infile)
labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = usage)
# print("labels: ", labels)
# ncol = len(datalist)
# print("ncol: ", ncol)
label_T, xT = datafile.FindDataArray(r'^T[\s|\(|\[|$]', flag = 'i')
label_sigma, ysigma = datafile.FindDataArray(r'^sigma[$|\S]?.*$', flag = 'i')
label_N, yN = datafile.FindDataArray(r'^N[$|\S]?.*$', flag = 'i')
ndata = len(xT)
print("ndata: ", ndata)
# print("xT=", label_T, xT)
# print("ysigma=", label_sigma, ysigma)
# print("yN=", label_N, yN)
print("Calculat EF, L and kappa,e")
yEF = []
yL = []
ykappa = []
print("kappa,e=sigma,obs*L,cal*T")
print("{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}"
.format("T(K)", "EF(eV)", "sigma(S/cm)", "Ne(cm^-3)", "L(Wohm/K^2)", "kappa,e"))
for i in range(len(xT)):
T = xT[i]
n = yN[i]
s = ysigma[i]
if n is not None:
EF, diffEF, flag = dos.EF_from_electrondensity(n, T, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100)
if not flag:
app.terminate("Error in calL(): EF calculation did not converge. diffEF={} eps={} in {} iteration".format(diffEF, epsEF, maxiter),
usage = usage, pause = True)
# print("n=", n, EF)
sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
= dos.cal_transport_properteis(T, EF, mobility, klatt, validate_error_str = 'Error in properteis()')
kappa = L * s / 0.01 * T
yEF.append(EF)
yL.append(L)
ykappa.append(kappa)
print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}"
.format(xT[i], EF, ysigma[i], yN[i], L, kappa))
header, ext = os.path.splitext(infile)
filebody = os.path.basename(header)
outxlsLfile = f'{filebody}-L-me{dos.meeff}-r{mobility.rfac}.xlsx'
print("")
print("Save data to [{}]".format(outxlsLfile))
tkVariousData().to_excel(outxlsLfile, ["T(K)", "EF,cal(eV)", "sigma,obs(S/cm)", "Ne,obs(cm^-3)", "L,cal(Wohm/K^2)", "kappa,e,cal=sigma,obs*L,cal*T(W/m/K)"],
[xT, yEF, ysigma, yN, yL, ykappa])
#=============================
# グラフの表示
#=============================
print("")
fig = plt.figure(figsize = figsize_sim)
if min(xT) == max(xT):
xX = yN
xlabel = '$N$ (cm$^{-3}$)'
xscale = 'log'
else:
xX = xT
xlabel = '$T$ (K)'
xscale = 'linear'
axEF = fig.add_subplot(2, 2, 1)
axNe = fig.add_subplot(2, 2, 2)
axsigma = axNe.twinx()
axL = fig.add_subplot(2, 2, 3)
axkappa = axL.twinx()
axkappa = fig.add_subplot(2, 2, 4)
axEF.plot(xX, yEF, label = 'EF(cal)', linestyle = '-', linewidth = 1.0, marker = 'o', markersize = 5.0)
axEF.set_xlabel(xlabel, fontsize = fontsize)
axEF.set_ylabel("$E_F$ (eV)", fontsize = fontsize)
axEF.set_xscale(xscale)
axEF.legend(fontsize = legend_fontsize, loc = 'best')
axEF.tick_params(labelsize = fontsize)
ins1 = axNe.plot (xX, yN, label = '$N_e$(obs)', linestyle = '-', linewidth = 1.0, color = 'red', marker = 'o', markersize = 5.0)
ins2 = axsigma.plot(xX, ysigma, label = r'$\sigma$(obs)', linestyle = 'dashed', linewidth = 1.0, color = 'blue', marker = 's', markersize = 5.0)
axNe.set_xlabel(xlabel, fontsize = fontsize)
axNe.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
axNe.set_xscale(xscale)
axNe.set_yscale('log')
axsigma.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
axsigma.set_yscale('log')
ins = ins1 + ins2
axNe.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
axNe.tick_params(labelsize = fontsize)
axsigma.tick_params(labelsize = fontsize)
axL.plot(xX, yL, label = 'L(cal)', linestyle = '-', linewidth = 1.0, marker = 'o', markersize = 5.0)
axL.set_xlabel(xlabel, fontsize = fontsize)
axL.set_ylabel(r"$L$ (W$\Omega$/K$^2$)", fontsize = fontsize)
axL.set_xscale(xscale)
axL.legend(fontsize = legend_fontsize, loc = 'best')
axL.tick_params(labelsize = fontsize)
axkappa.plot(xX, ykappa, label = r'$\kappa_e$(cal)', linestyle = '-', linewidth = 1.0, marker = 'o', markersize = 5.0)
axkappa.set_xlabel(xlabel, fontsize = fontsize)
axkappa.set_ylabel(r"$\kappa_e$=$\sigma$$L$$T$ (W/m/K)", fontsize = fontsize)
axkappa.set_xscale(xscale)
axkappa.set_yscale('log')
axkappa.legend(fontsize = legend_fontsize, loc = 'best')
axkappa.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
app.terminate("", pause = True)
[ドキュメント]
def calF():
"""
概要:
実測のホールキャリア濃度とホール移動度からホール因子を計算する。
詳細説明:
入力ファイルから温度、電気伝導度、ホールキャリア濃度、ホール移動度のデータを読み込み、
各データポイントにおけるフェルミ準位を計算する。
そのフェルミ準位と現在のモデルパラメータを用いてホール因子 `FH`、
ドリフトキャリア濃度 `n_drift`、ドリフト移動度 `mu_drift` を計算し、
結果をExcelファイル (`outxlsLfile`) に出力し、グラフ表示する。
引数:
なし。
"""
print("")
print("infile :", infile)
print("T_label : ", T_label)
print("n_label : ", n_label)
print("mu_label : ", mu_label)
print("sigma_label: ", sigma_label)
print("Carrier:")
print(" q={} e".format(mobility.charge))
print(" meff={}me".format(dos.meeff))
dos.NC = meff2NC_FEA(dos.meeff, T0)
dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
if mobility.charge < 0.0:
print(" NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
print(" DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
else:
print(" NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
print(" DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
print(" scattering factor (non-degenerated) r={}".format(mobility.rfac))
print(" mean free path prefactor: l0={} m".format(mobility.l0))
print("")
print("Read T, Ne, and sigma data from {}".format(infile))
datafile = tkVariousData(infile)
labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = usage)
label_T, xT = datafile.FindDataArray(T_label, flag = 'i')
label_sigma, ysigma = datafile.FindDataArray(sigma_label, flag = 'i')
label_N, yN = datafile.FindDataArray(n_label, flag = 'i')
label_mu, ymu = datafile.FindDataArray(mu_label, flag = 'i')
ndata = len(xT)
print("ndata: ", ndata)
print("")
print("Calculat EF, FH, mu_drift, and Ne")
yEF = []
yFH = []
yRH = []
yndrift = []
ymudrift = []
print("kappa,e=sigma,obs*L,cal*T")
print("{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}"
.format("T(K)", "EF(eV)", "sigma(S/cm)", "Ne,Hall(obs)", "mu,Hall(obs)", "Ne(cm-3)", "mu,drift(cm2/Vs)", "FH", "RH(cm3/C)"))
for i in range(len(xT)):
T = xT[i]
n = yN[i]
s = ysigma[i]
if n is not None:
EF, diffEF, flag = dos.EF_from_electrondensity(n, T, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100)
if not flag:
app.terminate("Error in calF(): EF calculation did not converge. diffEF={} eps={} in {} iteration".format(diffEF, epsEF, maxiter),
usage = usage, pause = True)
# sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
# = dos.cal_transport_properteis(T, EF, mobility, klatt, validate_error_str = 'Error in properteis()')
inf = dos.cal_Hall_properteis(T0, EF, mobility, validate_error_str = 'Error in properteis()')
FH = inf["FH"]
RH = inf["RH"]
RH0 = inf["RH0"]
n_drift = inf["nHall"]
mu_drift = inf["muHall"]
yEF.append(EF)
yFH.append(FH)
yRH.append(RH)
yndrift.append(yN[i] * FH)
ymudrift.append(ymu[i] / FH)
print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}"
.format(xT[i], EF, ysigma[i], yN[i], ymu[i], yndrift[i], ymudrift[i], FH, RH))
header, ext = os.path.splitext(infile)
filebody = os.path.basename(header)
outxlsLfile = f'{filebody}-FH-me{dos.meeff}-r{mobility.rfac}.xlsx'
print("")
print("Save data to [{}]".format(outxlsLfile))
tkVariousData().to_excel(outxlsLfile, ["T(K)", "EF,cal(eV)", "sigma,obs(S/cm)", "Ne,Hall(obs)(cm-3)", "mu,Hall(obs)(cm2/Vs)",
"Ne=Ne,Hall*FH(cm-3)", "mu,drift=mu,Hall/FH(cm2/Vs)", "FH", "RH(cm3/C)"],
[xT, yEF, ysigma, yN, ymu, yndrift, ymudrift, yFH, yRH])
#=============================
# グラフの表示
#=============================
print("")
fig = plt.figure(figsize = figsize_sim)
if min(xT) == max(xT):
xX = yN
xlabel = '$N$ (cm$^{-3}$)'
xscale = 'log'
else:
xX = xT
xlabel = '$T$ (K)'
xscale = 'linear'
axEF = fig.add_subplot(2, 3, 1)
axsigma = fig.add_subplot(2, 3, 2)
axNe = fig.add_subplot(2, 3, 3)
axmu = fig.add_subplot(2, 3, 4)
axFH = fig.add_subplot(2, 3, 5)
axRH = fig.add_subplot(2, 3, 6)
axEF.plot(xX, yEF, label = 'EF(cal)', linestyle = '-', linewidth = 1.0, marker = 'o', markersize = 5.0)
axEF.set_xlabel(xlabel, fontsize = fontsize)
axEF.set_ylabel("$E_F$ (eV)", fontsize = fontsize)
axEF.set_xscale(xscale)
axEF.legend(fontsize = legend_fontsize, loc = 'best')
axEF.tick_params(labelsize = fontsize)
axsigma.plot(xX, ysigma, label = r'$\sigma$(obs)', linestyle = 'dashed', linewidth = 1.0, color = 'blue', marker = 's', markersize = 5.0)
axsigma.set_xlabel(xlabel, fontsize = fontsize)
axsigma.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
axsigma.set_xscale(xscale)
axsigma.set_yscale('log')
axsigma.legend(fontsize = legend_fontsize, loc = 'best')
axsigma.tick_params(labelsize = fontsize)
axNe.plot (xX, yN, label = '$N_e$$_{,Hall}(obs)', linestyle = '-', linewidth = 1.0, color = 'red', marker = 'o', markersize = 5.0)
axNe.plot (xX, yndrift, label = '$N_e$$_{,drift}$(obs)', linestyle = '-', linewidth = 0.5, color = 'green', marker = '^', markersize = 3.0)
axNe.set_xlabel(xlabel, fontsize = fontsize)
axNe.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
axNe.set_xscale(xscale)
axNe.set_yscale('log')
axNe.legend(fontsize = legend_fontsize, loc = 'best')
axNe.tick_params(labelsize = fontsize)
axmu.plot(xX, ymu, label = r'$\mu_{Hall}$(obs)', linestyle = '-', linewidth = 1.0, color = 'red', marker = 'o', markersize = 5.0)
axmu.plot(xX, ymudrift, label = r'$\mu_{drift}$(obs)', linestyle = '-', linewidth = 0.5, color = 'green', marker = '^', markersize = 3.0)
axmu.set_xlabel(xlabel, fontsize = fontsize)
axmu.set_ylabel(r"$\mu$ (cm$^2$/Vs)", fontsize = fontsize)
axmu.set_xscale(xscale)
axmu.legend(fontsize = legend_fontsize, loc = 'best')
axmu.tick_params(labelsize = fontsize)
axFH.plot(xX, yFH, label = '$F_{Hall}$', linestyle = '-', linewidth = 1.0, color = 'blue', marker = 'o', markerfacecolor = 'blue', markersize = 3.0)
axFH.set_xlabel(xlabel, fontsize = fontsize)
axFH.set_ylabel("$F_{Hall}$", fontsize = fontsize)
axFH.set_xscale(xscale)
axFH.legend(fontsize = legend_fontsize, loc = 'best')
axFH.tick_params(labelsize = fontsize)
axRH.plot(xX, yRH, label = '$R_H$', linestyle = '-', linewidth = 1.0, color = 'blue', marker = 'o', markerfacecolor = 'blue', markersize = 3.0)
axRH.set_xlabel(xlabel, fontsize = fontsize)
axRH.set_ylabel("$R_H$ (cm$^3$/C)", fontsize = fontsize)
axRH.set_xscale(xscale)
axRH.set_yscale('log')
axRH.legend(fontsize = legend_fontsize, loc = 'best')
axRH.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
app.terminate("", pause = True)
[ドキュメント]
def T():
"""
概要:
温度依存性を考慮した熱電特性のシミュレーションとプロットを行う。
詳細説明:
入力ファイルから実測データを読み込み、指定された温度範囲 `Tmin` から `Tmax` で、
キャリア濃度 `N` に対する様々な熱電特性（フェルミ準位 `EF`、ゼーベック係数 `S`、
電気伝導度 `sigma`、移動度 `mu`、電子熱伝導率 `kappa`、出力因子 `PF`、ZT値 `ZT`）の
計算を行う。結果はExcelファイルに保存され、温度をパラメータとしたプロットとして表示される。
引数:
なし。
"""
global label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu
print("")
print("Carrier:")
print(" q={} e".format(mobility.charge))
print(" meff={}me".format(dos.meeff))
dos.NC = meff2NC_FEA(dos.meeff, T0)
dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
if mobility.charge < 0.0:
print(" NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
print(" DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
else:
print(" NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
print(" DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
print(" scattering factor (non-degenerated) r={}".format(mobility.rfac))
print(" mean free path prefactor: l0={} m".format(mobility.l0))
print("")
print("Read S data from {}".format(infile))
label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu \
= read_datafile(infile, usage = usage)
label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu \
= construct_lists(label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu)
print("x=", xsample)
ndata = len(xsample)
print("ndata: ", ndata)
if yN[0] is not None:
print("{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:10}\t{:>10}"
.format("sample", "EF(eV)", "S(uV/K)", "sigma(S/cm)", "mu(cm2/Vs)", "Ne(cm-3)"))
else:
print("{:>10}\t{:>10}\t{:10}\t{:>10}".format("sample", "EF(eV)", "S(uV/K)", "sigma(S/cm)"))
for i in range(ndata):
sample = xsample[i]
S = yS[i]
sigma = ysigma[i]
n = yN[i]
mu = ymu[i]
if n is not None:
print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, S, sigma, mu, n))
else:
print("{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, S, sigma))
Nmin_data = min(yN)
Nmax_data = max(yN)
if Nmin_data < Nmin:
Nmin = Nmin_data
if Nmax < Nmax_data:
Nmax = Nmax_data
print("")
print(f"Plot N range: {Nmin:12.4g} - {Nmax:12.4g}")
lnNmin = log(Nmin)
lnNmax = log(Nmax)
lnNstep = (lnNmax - lnNmin) / (nN - 1)
print("")
print("N range: {:10.4g} - {:10.4g} cm-3, {:10.4g} step in ln(N), nN={}".format(Nmin, Nmax, lnNstep, nN))
xEFcal = make_matrix2(nT, nN, 0.0)
xNsim = make_matrix2(nT, nN, 0.0)
ysigmasim = make_matrix2(nT, nN, 0.0)
ymusim = make_matrix2(nT, nN, 0.0)
yScal = make_matrix2(nT, nN, 0.0)
ykappacal = make_matrix2(nT, nN, 0.0)
yPFcal = make_matrix2(nT, nN, 0.0)
yZTcal = make_matrix2(nT, nN, 0.0)
Tstep = (Tmax - Tmin) / (nT - 1)
for iT in range(nT):
T = Tmin + iT * Tstep
dos.NC = meff2NC_FEA(dos.meeff, T)
dos.DC0 = meff2DC0_FEA(dos.meeff, T)
print("{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}"
.format("T(K)", "N(cm-3)", "EF(eV)", "S,cal(uV/K)", "sigma,cal(S/cm)", "Ne,cal(cm-3)", "mu,cal(cm2/Vs)", "kappa,e,cal(W/m/K)", "PF(uW/cm/K^2)", "ZT"))
for iN in range(nN):
lnN = lnNmin + lnNstep * iN
N = exp(lnN)
# print("N=", N, dos.NC)
if N < dos.NC:
EF = kB * T0 * log(N / dos.NC) / e
else:
EF = BMShift_FEA(N, dos.DC0)
EF, diffEF, flag = dos.EF_from_electrondensity(N, T, EF0 = EF, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100)
if not flag:
app.terminate("Error in sim(): EF calculation did not converge. diffEF={} eps={} in {} iteration".format(diffEF, epsEF, maxiter),
usage = usage, pause = True)
sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
= dos.cal_transport_properteis(T, EF, mobility, klatt, validate_error_str = 'Error in sim()')
xEFcal[iT][iN] = EF
xNsim[iT][iN] = n
ysigmasim[iT][iN] = sigma
ymusim[iT][iN] = mu
yScal[iT][iN] = S
ykappacal[iT][iN] = kappa
yPFcal[iT][iN] = PF
yZTcal[iT][iN] = ZT
print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}"
.format(T, N, EF, S, sigma, n, mu, kappa, PF, ZT))
print("")
print("Save parameters to [{}]".format(parameterfile))
save_parameterfile(S2 = '')
print("")
labels = ["N(cm-3)"] + ["{:.4g} K".format(Tmin + iT * Tstep) for iT in range(nT)]
header, ext = os.path.splitext(infile)
filebody = os.path.basename(header)
xlsfile = filebody + '-EF.xlsx'.format(T)
print("Save data to [{}]".format(xlsfile))
# print("labels=", labels)
# print("datalist=", [xNsim[0]])# + [xEFcal[iT] for iT in range(nT)])
tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + xEFcal) #[xEFcal[iT] for iT in range(nT)])
xlsfile = filebody + '-sigma.xlsx'.format(T)
print("Save data to [{}]".format(xlsfile))
tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + ysigmasim)
xlsfile = filebody + '-mu.xlsx'.format(T)
print("Save data to [{}]".format(xlsfile))
tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + ymusim)
xlsfile = filebody + '-S.xlsx'.format(T)
print("Save data to [{}]".format(xlsfile))
tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + yScal)
xlsfile = filebody + '-kappa.xlsx'.format(T)
print("Save data to [{}]".format(xlsfile))
tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + ykappacal)
xlsfile = filebody + '-PF.xlsx'.format(T)
print("Save data to [{}]".format(xlsfile))
tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + yPFcal)
xlsfile = filebody + '-ZT.xlsx'.format(T)
print("Save data to [{}]".format(xlsfile))
tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + yZTcal)
#=============================
# グラフの表示
#=============================
print("")
fig = plt.figure(figsize = figsize)
axSsobs = fig.add_subplot(3, 3, 1)
axSsobsb = axSsobs.twinx()
axmuNobs = fig.add_subplot(3, 3, 2)
axmuNobsb = axmuNobs.twinx()
axNs = fig.add_subplot(3, 3, 3)
axNmu = fig.add_subplot(3, 3, 4)
axNS = fig.add_subplot(3, 3, 5)
axNkappa = fig.add_subplot(3, 3, 6)
axNPF = fig.add_subplot(3, 3, 7)
axNZT = fig.add_subplot(3, 3, 8)
axNEF = fig.add_subplot(3, 3, 9)
maxS = max(yS)
minS = min(yS)
if maxS > 0.0:
maxS *= 1.2
else:
maxS = 0.0
if minS > 0.0:
minS = 0.0
else:
minS *= 1.2
ins1 = axSsobs.plot (xsample, yS, label = 'S', linestyle = '-', linewidth = 1.0, color = 'red', marker = 'o', markersize = 5.0)
ins2 = axSsobsb.plot(xsample, ysigma, label = r'$\sigma$', linestyle = '-', linewidth = 1.0, color = 'blue', marker = '^', markersize = 10.0)
axSsobs.set_xlabel("sample", fontsize = fontsize)
axSsobs.set_ylabel(r"S (K/$\mu$V)", fontsize = fontsize)
axSsobsb.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
ins = ins1 + ins2
axSsobs.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
axSsobs.tick_params(labelsize = fontsize)
axSsobsb.tick_params(labelsize = fontsize)
if yN[0] is not None:
ins1 = axmuNobs.plot( xsample, ymu, label = r'$\mu$', linestyle = '-', color = 'red', linewidth = 0.5, marker = 'o', markersize = 10.0)
ins2 = axmuNobsb.plot(xsample, yN, label = '$N_e$', linestyle = '-', color = 'blue', linewidth = 0.5, marker = 'o', markersize = 10.0)
axmuNobs.set_xlabel("sample", fontsize = fontsize)
axmuNobs.set_ylabel(r"$\mu$ (cm$^2$/V/s)", fontsize = fontsize)
axmuNobsb.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
axmuNobsb.set_yscale('log')
ins = ins1 + ins2
axmuNobs.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
axmuNobs.tick_params(labelsize = fontsize)
axmuNobsb.tick_params(labelsize = fontsize)
if yN[0] is not None:
axNs.plot( yN, ysigma, label = r'$\sigma_{,obs}$', linestyle = '', marker = 'o', markersize = 5.0)
axNmu.plot(yN, ymu, label = '$N_e$$_{,obs}$', linestyle = '', marker = 'o', markersize = 5.0)
axNS.plot( yN, yS, label = '$S_{obs}$', linestyle = '', marker = 'o', markersize = 5.0)
for iT in range(nT):
T= Tmin + iT * Tstep
axNs.plot( xNsim[iT], ysigmasim[iT], label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
axNmu.plot( xNsim[iT], ymusim[iT], label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
axNS.plot( xNsim[iT], yScal[iT], label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
axNkappa.plot(xNsim[iT], ykappacal[iT], label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
axNPF.plot( xNsim[iT], yPFcal[iT], label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
axNZT.plot( xNsim[iT], yZTcal[iT], label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
axNEF.plot( xNsim[iT], xEFcal[iT], label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
axNs.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
axNmu.set_xlabel(r"$N$ (cm$^{-3}$)", fontsize = fontsize)
axNS.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
axNkappa.set_xlabel( "$N$ (cm$^{-3}$)", fontsize = fontsize)
axNPF.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
axNZT.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
axNEF.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
axNs.set_ylabel( r"$\sigma$ (S/cm)", fontsize = fontsize)
axNmu.set_ylabel(r"$\mu$ (cm$^2$/V/s)", fontsize = fontsize)
axNS.set_ylabel( "S (V/K)", fontsize = fontsize)
axNkappa.set_ylabel( r"$\kappa$$_{,e}$ (W/m/K)", fontsize = fontsize)
axNPF.set_ylabel( r"PF ($\mu$W/cm/K$^2$)", fontsize = fontsize)
axNZT.set_ylabel( "ZT", fontsize = fontsize)
axNEF.set_ylabel( "$E_F$", fontsize = fontsize)
axNs.set_xscale('log')
axNmu.set_xscale('log')
axNS.set_xscale('log')
axNkappa.set_xscale('log')
axNPF.set_xscale('log')
axNZT.set_xscale('log')
axNEF.set_xscale('log')
axNs.set_yscale('log')
axNkappa.set_yscale('log')
axNs.legend(fontsize = legend_fontsize, loc = 'best')
axNmu.legend(fontsize = legend_fontsize, loc = 'best')
axNS.legend(fontsize = legend_fontsize, loc = 'best')
axNkappa.legend(fontsize = legend_fontsize, loc = 'best')
axNPF.legend(fontsize = legend_fontsize, loc = 'best')
axNZT.legend(fontsize = legend_fontsize, loc = 'best')
axNEF.legend(fontsize = legend_fontsize, loc = 'best')
axNs.tick_params(labelsize = fontsize)
axNmu.tick_params(labelsize = fontsize)
axNS.tick_params(labelsize = fontsize)
axNkappa.tick_params(labelsize = fontsize)
axNPF.tick_params(labelsize = fontsize)
axNZT.tick_params(labelsize = fontsize)
axNEF.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
app.terminate("", pause = True)
[ドキュメント]
def main():
"""
概要:
スクリプトの主要な実行フローを管理する。
詳細説明:
アプリケーションの初期化、ログファイルのリダイレクト、コマンドライン引数の解析と
グローバル変数の更新、INIファイルからのパラメータ読み込み、
そして指定された `mode` に応じた各機能関数の呼び出しを行う。
処理の開始時と終了時に情報をコンソールに出力する。
引数:
なし。
"""
global app
global parameterfile, parameterbkfile, outxlsfile, outxlsPropfile
global outxlsfFDfile, outxlsFjfile
app = tkApplication()
logfile = app.replace_path(infile)
print(f"Open logfile [{logfile}]")
app.redirect(targets = ["stdout", logfile], mode = 'w')
updatevars()
parameterfile = app.replace_path(infile, template = ["{dirname}", "{filebody}.in"])
parameterbkfile = app.replace_path(infile, template = ["{dirname}", "{filebody}.in.bak"])
outxlsfile = app.replace_path(infile, template = ["{dirname}", "{filebody}-Seebeck-{mode}.xlsx"], ext_dict = {"mode": mode})
outxlsPropfile = app.replace_path(infile, template = ["{dirname}", "{filebody}-properties-{mode}.xlsx"], ext_dict = {"mode": mode})
outxlsfFDfile = app.replace_path(infile, template = ["{dirname}", "fFD.xlsx"])
outxlsFjfile = app.replace_path(infile, template = ["{dirname}", "Fj.xlsx"])
print("")
print("mode: ", mode)
print("infile : {}".format(infile))
print("out xlsx file : {}".format(outxlsfile))
print("parameter file : {}".format(parameterfile))
print("parameter backup file: {}".format(parameterbkfile))
print("")
print("")
print("Read [{}]".format(parameterfile))
read_parameters(parameterfile)
# check args again to use given parameters
updatevars()
print("")
print("Electronic structure")
# print(f"Eg={dos.EC - dos.EV} eV")
# print(f"EV={dos.EV} eV")
print(f"EC={dos.EC} eV")
# print(f"EA={dos.EA} eV")
# print(f"NA={dos.NA} cm^-3")
# print(f"ED={dos.ED} eV")
# print(f"ND={dos.ND} cm^-3")
# print(f"mh_eff={dos.mheff} me_0")
print(f"me_eff={dos.meeff} me_0")
print(f"Initial parameter")
print(f" EF0={dos.EF0} eV")
print("Transport parameters")
print(f"charge={mobility.charge} e")
print(f"meff={mobility.meff} me_0")
print(f"rfac={mobility.rfac}")
print(f"l0 ={mobility.l0} m")
if mode == 'basic':
basic()
elif mode == 'prop':
properties()
elif mode == 'Hall':
Hall()
elif mode == 'init':
init()
elif mode == 'sim':
sim()
elif mode == 'fit':
fit()
elif mode == 'T':
T()
elif mode == 'calL':
calL()
elif mode == 'calF':
calF()
else:
app.terminate("Error in main: Invalide mode [{}]".format(mode), usage = usage, pause = True)
if __name__ == "__main__":
main()