"""
N-T Fit Semi-FEA: ホール効果データを用いた半導体キャリア濃度・移動度解析スクリプト
概要:
温度依存性のあるホール効果データ（キャリア濃度、移動度、ホール係数など）を
自由電子近似 (FEA) モデルに基づいて解析し、バンド構造や不純物準位に関する
パラメータをフィッティング、またはシミュレーションします。
詳細説明:
このスクリプトは、測定されたホール効果の温度依存性データ（ホールキャリア濃度、
ホール移動度、導電率、ホール係数など）を読み込み、tkDOS (半導体統計計算) および
tkMobility (移動度計算) モジュールと連携して、自由電子近似モデルを用いた
理論曲線と実験データを比較・解析します。
以下のモードをサポートします。
- 'sim': 現在のパラメータでキャリア濃度、移動度、ホール係数などの温度依存性を
シミュレーションし、結果をプロットして保存します。
- 'fit' / 'fit_e' / 'fit_h': ホールキャリア濃度データ (`Ns = 1/eRH`) に対して
非線形最小二乗フィッティングを行います。
'fit' は二キャリアモデル、'fit_e' は電子のみ、
'fit_h' はホールのみの単一キャリアモデルを使用します。
散乱因子は `r=0.5` に固定されます。
- 'fit_RH': ホール係数 (`RH`) データに対して非線形最小二乗フィッティングを行います。
二キャリアモデルを使用し、散乱因子 `r` は `tkMobility` オブジェクトで
設定された値を使用します。
最適化には `scipy.optimize.minimize` を使用し、指定されたアルゴリズムで
DOSパラメータ（有効質量、バンドギャップ、不純物準位、不純物濃度など）を調整します。
フィッティングの進行状況はリアルタイムでコンソール出力およびグラフ更新によって
可視化されます。
関連リンク:
:doc:`N-T-Fit-semi_FEA_usage`
"""
import os
import sys
import signal
import numpy as np
from numpy import sin, cos, tan, pi, exp, log
from matplotlib import pyplot as plt
from scipy import optimize # newton関数はscipy.optimizeモジュールに入っている
from scipy.optimize import minimize
sys.path.append("c:/tkProg/tklib/python")
sys.path.append("d:/tkProg/tklib/python")
from tklib.tkfile import tkFile
from tklib.tkutils import terminate, getarg, getintarg, getfloatarg, save_csv
from tklib.tkutils import is_exist, is_file, is_dir
#from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg
from tklib.tkapplication import tkApplication
from tklib.tksci.tksci import pi, h, hbar,c, e, kB, me, Ea_Arrhenius
from tklib.tktransport.tkTransport import read_Hall_excel
from tklib.tktransport.tkTransport import integrate_Simpson_list, fe, fh, NC2meff_FEA
from tklib.tktransport.tkDOS_FEA import integrate_Simpson, integrate_Simpson_list, tkDOS
from tklib.tktransport.tkmobility_tau import tkMobility, split_optstr
usage_str = '''
""
"Usage: Variables in () are optional"
" python {} mode (file args)".format(sys.argv[0])
" (i) python {} sim infile(.xlsx) EC EA NA ED ND".format(sys.argv[0])
" Plot N-T with simulated data"
" ex: python {} sim {} {} {} {} {} {}".format(sys.argv[0], cparams.infile, dos.EC, dos.EA, dos.NA, dos.ED, dos.ND)
" (ii) python {} fit|fit_h|fit_e infile(.xlsx) EC EA NA ED ND Eb a1 a2)".format(sys.argv[0])
" Fit to N-mu-T data"
" mode = fit : Fit to Ns = 1/eRH by two-carrier (hole-electron) model. Scattering factor is fixed to r = 0.5"
" fit_h : Fit to Nh by single-carrier model. Scattering factor is fixed to r = 0.5"
" fit_e : Fit to Ne by single-carrier model. Scattering factor is fixed to r = 0.5"
" fit_RH: Fit to RH by two-carrier (hole-electron) model with scattering factor"
" ex: python {} fit {} {} {} {} {} {}".format(sys.argv[0], cparams.infile, dos.EC, dos.EA, dos.NA, dos.ED, dos.ND)
'''[1:-1]
[ドキュメント]
def pint(str):
"""
概要: 文字列を整数に変換します。
詳細説明:
Python組み込みのint()関数を呼び出します。
この関数はtklib.tkutilsにも同名の関数がありますが、ここではローカルで再定義されています。
:param str: 整数に変換する文字列。
:type str: str
:returns: 変換された整数値。
:rtype: int
"""
return int(str)
[ドキュメント]
def pfloat(str):
"""
概要: 文字列を浮動小数点数に変換します。
詳細説明:
Python組み込みのfloat()関数を呼び出します。
この関数はtklib.tkutilsにも同名の関数がありますが、ここではローカルで再定義されています。
:param str: 浮動小数点数に変換する文字列。
:type str: str
:returns: 変換された浮動小数点数値。
:rtype: float
"""
return float(str)
[ドキュメント]
def sigint_handler(signal, frame):
"""
概要: Ctrl+C (SIGINT) シグナルを処理するハンドラ関数です。
詳細説明:
キーボードからの割り込み (Ctrl+C) を検知した際に呼び出されます。
"KeyboardInterrupt is caught" というメッセージを表示し、
アプリケーションを強制終了します。
:param signal: シグナル番号 (SIGINT)。
:type signal: int
:param frame: 現在のスタックフレーム。
:type frame: types.FrameType
:returns: None
:rtype: None
"""
print ('KeyboardInterrupt is caught')
terminate(usage = usage)
signal.signal(signal.SIGINT, sigint_handler)
# EPS to avoid log(0.0) => log(x + cparams.eps)
eps = 1.0e-700
[ドキュメント]
def initialize():
"""
概要: アプリケーションの初期設定を行い、主要なオブジェクトを初期化します。
詳細説明:
`tkApplication`, `tkDOS`, `tkMobility` の各オブジェクトを生成し、
それらのデフォルトパラメータを設定します。これには、フィッティングモード、
入力ファイルパス、半導体のDOS (状態密度) パラメータ（有効質量、バンドギャップ、
不純物準位と濃度）、および移動度パラメータ（散乱因子など）が含まれます。
最適化アルゴリズムやグラフ表示に関する設定も行われます。
:returns: 初期化されたアプリケーション、パラメータ、DOS、移動度の各オブジェクト。
:rtype: tuple[tkApplication, tkParams, tkDOS, tkMobility]
"""
#================================
# Global variables
#================================
app = tkApplication(usage_str = usage_str)
argv, narg = app.get_argv()
app.dos = tkDOS()
app.mu = tkMobility()
cparams = app.get_params()
cparams.debug = 0
cparams.use_simple = 0
cparams.print_level = 0
# mode: 'sim', 'fit', 'fit_e', 'fit_h', "fit_RH"
cparams.mode = 'sim'
# cparams.mode = 'fit'
# cparams.mode = 'fit_e'
# cparams.mode = 'fit_h'
# data file to fit (.xlsx)
cparams.infile = None
cparams.parameterfile = None
cparams.parameterbkfile = None
cparams.outcsvfile = None
dos = app.dos
# Assume me = mh to calculate RH for two-carrier model
dos.cal_mobility_e = lambda T: 1.0
dos.cal_mobility_h = lambda T: 1.0
dos.cal_mobility_e_description = "Same mobility for hole and electron"
dos.cal_mobility_h_description = "Same mobility for hole and electron"
dos.meeff = 0.3
dos.mheff = 1.0
dos.EV = 0.0
dos.EC = 1.1
dos.EA = 0.05 # Acceptor level, eV
dos.NA = 0.0e17
dos.ED = 1.05 # Donor level
dos.ND = 0.0e17
cparams.EF0 = 0.0 # initial EF to find EF
#移動度パラメータ
mobility = app.mu
# For debug purpose. 1 will use 3-terms polynomial for the inverse of mobility
mobility.debug = cparams.debug
mobility.use_simple = cparams.use_simple
#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
dos.varkeys = ["meeff", "mheff", "EV", "EC", "EA", "NA", "ED", "ND"]
# For mode = 'me'. E range for fitting to me*(DOS)
cparams.dE_me_lsq = 0.5 # measured from EC/EV, eV
# EF range for mode = 'EF'
cparams.T0 = 300 # K
cparams.dEFmin = -1.0 # measured from EV, eV
cparams.dEFmax = 1.0 # measured from EC, eV
cparams.nEF = 50
cparams.Estep = 0.01 # Integration step, eV
# Temperature range for mode = 'T'
cparams.Tmin = 300 # K
cparams.Tmax = 600
cparams.nT = 11
cparams.Tstep = None
# Integration range w.r.t. kBT
cparams.Einteg0 = -3.0
cparams.Einteg1 = 3.0
cparams.nrange = 6.0
# Bisection method parameters
# Initial search range
cparams.dE_bisec = 0.4
# EPS, max iteration number, Output cycle
cparams.eps_bisec = 1.0e-7
cparams.nmaxiter_bisec = 200
cparams.iprintiterval_bisec = 1
# 最適化パラメータの初期値
dos.varname = ["meeff", "mheff", "EV", "EC", "EA", "NA", "ED", "ND"]
dos.varunit = [ "me", "me", "eV", "eV", "eV", "cm-3", "eV", "cm-3"]
dos.ai0 = dos.parameter_list()
dos.optid = [ 1, 1, 0, 0, 1, 1, 0, 0]
dos.xT = []
dos.ys = []
dos.ymu = []
dos.yNhini = []
dos.yNeini = []
dos.yNsini = []
dos.yNsabsini = []
#=============================================
# 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)
cparams.method = "nelder-mead"
#method = 'cg'
#method = 'powell'
#method = 'bfgs'
cparams.maxiter = 300
cparams.tol = 1.0e-4
cparams.h_diff = 1.0e-3
cparams.outputinterval = 1
#=============================
# Graph configuration
#=============================
cparams.fig = None
cparams.figsize = [12, 8]
cparams.fontsize = 14
cparams.legend_fontsize = 12
cparams.graphupdateinterval = 10
return app, cparams, app.dos, app.mu
[ドキュメント]
def print_parameters(varname, varunit, ai, optid):
"""
概要: 指定されたパラメータリストを整形して標準出力に表示します。
詳細説明:
パラメータ名、その値、単位、およびそのパラメータが最適化の対象であるか
を示すID (`optid`) を含めて表示します。
:param varname: パラメータ名のリスト。
:type varname: list[str]
:param varunit: パラメータ単位のリスト。
:type varunit: list[str]
:param ai: パラメータ値のリスト。
:type ai: list[float]
:param optid: 各パラメータが最適化対象かを示すIDのリスト (1:対象, 0:非対象)。
:type optid: list[int]
:returns: None
:rtype: None
"""
for i in range(len(varname)):
print(" {:10s}: {:14.8g} {:10}: optid={}".format(varname[i], ai[i], varunit[i], optid[i]))
[ドキュメント]
def usage(app):
"""
概要: アプリケーションの使用方法メッセージを表示します。
詳細説明:
`tkApplication` オブジェクトに格納されている `usage_str` を解析し、
整形して標準出力に表示します。
:param app: アプリケーションオブジェクト。
:type app: tkApplication
:returns: None
:rtype: None
"""
# app.usage(infile = cparams.infile)
for s in app.usage_str.split('\n'):
cmd = 'print({})'.format(s.rstrip())
eval(cmd)
[ドキュメント]
def updatevars(app, cparams, dos, mu):
"""
概要: コマンドライン引数に基づいてアプリケーションの設定を更新します。
詳細説明:
sys.argvからフィッティングモード、入力ファイル名、DOSパラメータ（EC, EA, NA, ED, ND）を抽出し、
`cparams` と `dos` オブジェクトの対応する属性を更新します。
入力ファイルが指定されていない場合はエラーで終了します。
ファイルパスの拡張子処理も行います。
:param app: アプリケーションオブジェクト。
:type app: tkApplication
:param cparams: アプリケーションパラメータオブジェクト。
:type cparams: tkParams
:param dos: 状態密度オブジェクト。
:type dos: tkDOS
:param mu: 移動度オブジェクト (この関数内では直接使用されない)。
:type mu: tkMobility
:returns: None
:rtype: None
"""
argv = sys.argv
if len(argv) == 1:
usage(app)
exit()
cparams.mode = getarg( 1, cparams.mode)
cparams.infile = getarg( 2, cparams.infile)
dos.EC = getfloatarg(3, dos.EC)
dos.EA = getfloatarg(4, dos.EA)
dos.NA = getfloatarg(5, dos.NA)
dos.ED = getfloatarg(6, dos.ED)
dos.ND = getfloatarg(7, dos.ND)
dos.Eg = dos.EC - dos.EV
ai0 = dos.parameter_list()
if cparams.infile is None:
terminate("Error in updatevars: infile must be specified", usage = usage)
header, ext = os.path.splitext(cparams.infile)
filebody = os.path.basename(header)
cparams.infile = filebody + '.xlsx'
if not is_file(cparams.infile):
cparams.infile = filebody + '.csv'
cparams.parameterfile = filebody + '.in'
cparams.parameterbkfile = filebody + '.in.bak'
cparams.outcsvfile = filebody + f'-EFT-{cparams.mode}.csv'
[ドキュメント]
def recover_parameters(ais, optid, aidef0):
"""
概要: 最適化された可変パラメータと固定パラメータを結合し、完全なパラメータリストを再構築します。
詳細説明:
`optid` リストに基づいて、`ais` (最適化されたパラメータのサブセット) と
`aidef0` (最適化されない固定パラメータ、またはデフォルト値) を結合し、
元のパラメータリストと同じ構造で完全なリストを生成して返します。
:param ais: 最適化によって更新されたパラメータ値のリスト。
:type ais: list[float]
:param optid: 各パラメータが最適化対象かを示すIDのリスト (1:対象, 0:非対象)。
:type optid: list[int]
:param aidef0: 全パラメータのデフォルト値または固定値のリスト。
:type aidef0: list[float]
:returns: 完全なパラメータリスト。
:rtype: list[float]
"""
aidef = 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` リストで `1` とマークされたパラメータ (`ais0` 内の) のみを抽出し、
最適化関数 (`scipy.optimize.minimize`) に渡すための新しいリストを生成します。
:param ais0: 全パラメータ値のリスト。
:type ais0: list[float]
:param optid: 各パラメータが最適化対象かを示すIDのリスト (1:対象, 0:非対象)。
:type optid: list[int]
:returns: 最適化対象のパラメータのみを含むリスト。
:rtype: list[float]
"""
ais = ais0.copy()
#ai0 = [mu0, Eb, a1, a2]
ai = []
c = 0
for i in range(len(optid)):
if optid[i] == 1:
ai.append(ais[i])
return ai
[ドキュメント]
def calS2(app, cparams, dos, mobility, ai):
"""
概要: キャリア濃度データに対する二乗誤差 (S2) を計算します。
詳細説明:
入力されたパラメータ `ai` を使用して `tkDOS` オブジェクトのパラメータを更新し、
指定された温度リスト `dos.xT` における理論的なキャリア濃度 (`yNc`) を計算します。
観測されたキャリア濃度 `dos.yN` と計算された `yNc` の対数間の二乗誤差を計算し、
その平均値を返します。フィッティングモードに応じて、電子、ホール、または
絶対キャリア濃度が比較対象となります。
:param app: アプリケーションオブジェクト。
:type app: tkApplication
:param cparams: アプリケーションパラメータオブジェクト。
:type cparams: tkParams
:param dos: 状態密度オブジェクト。
:type dos: tkDOS
:param mobility: 移動度オブジェクト (この関数内では直接使用されない)。
:type mobility: tkMobility
:param ai: 現在の最適化パラメータ値のリスト。
:type ai: list[float]
:returns: 計算された二乗誤差 (S2) の値。
:rtype: float
"""
aiorg = dos.parameter_list()
ai0 = dos.parameter_list()
ainew = recover_parameters(ai, dos.optid, ai0)
dos.set_parameters(ainew)
xTinv, yEFfin, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, yNsabsfin, ydQ = dos.cal_T_lists(dos.xT, print_level = cparams.print_level)
# print("")
if cparams.print_level >= 1:
print("{:8}\t{:8}\t{:10}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}"
.format('T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,fin', 'Nh,fin', 'NA-,fin', 'ND+,fin', 'RH,fin', '|Ns,fin|'))
for i in range(len(dos.xT)):
print("{:8.4f}\t{:8.4f}\t{:10.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}"
.format(dos.xT[i], xTinv[i], yEFfin[i], dos.yN[i], yNefin[i], yNhfin[i], yNAmfin[i], yNDpfin[i], yRHfin[i], yNsabsfin[i]))
if cparams.mode == 'fit_h':
yNc = yNhfin
elif cparams.mode == 'fit_e':
yNc = yNefin
else:
yNc = yNsabsfin
S2 = 0.0
for iT in range(len(dos.xT)):
T = dos.xT[iT]
S2 += (log(dos.yN[iT]+eps) - log(yNc[iT]+eps))**2
S2 /= (len(dos.xT) - 1)
dos.set_parameters(aiorg)
return S2
[ドキュメント]
def calS2RH(app, cparams, dos, mobility, ai):
"""
概要: ホール係数 (RH) またはキャリア濃度データに対する二乗誤差 (S2) を計算します。
詳細説明:
入力されたパラメータ `ai` を使用して `tkDOS` オブジェクトのパラメータを更新し、
指定された温度リスト `dos.xT` における理論的なキャリア濃度やホール係数を計算します。
フィッティングモードが `fit_RH` の場合は、観測されたホール係数 `dos.yRH` と
計算された `yRHfin` の対数間の二乗誤差を計算します。それ以外のモードでは、
`calS2` と同様にキャリア濃度間の対数誤差を計算します。
:param app: アプリケーションオブジェクト。
:type app: tkApplication
:param cparams: アプリケーションパラメータオブジェクト。
:type cparams: tkParams
:param dos: 状態密度オブジェクト。
:type dos: tkDOS
:param mobility: 移動度オブジェクト。
:type mobility: tkMobility
:param ai: 現在の最適化パラメータ値のリスト。
:type ai: list[float]
:returns: 計算された二乗誤差 (S2) の値。
:rtype: float
"""
aiorg = dos.parameter_list()
ai0 = dos.parameter_list()
ainew = recover_parameters(ai, dos.optid, ai0)
dos.set_parameters(ainew)
# xTinv, yEFfin, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, yNsabsfin, ydQ = dos.cal_T_lists(dos.xT, print_level = cparams.print_level)
inf = dos.cal_T_lists2(dos.xT, mobility, print_level = cparams.print_level)
xTinv = inf['Tinv']
yEFfin = inf['EF']
yNefin = inf['n_e']
yNhfin = inf['n_h']
# yNeHallfin = inf['nHall_e']
# yNhHallfin = inf['nHall_h']
yNsfin = inf['n_s']
yNsabsfin = inf['n_s(abs)']
yNAmfin = inf['NAm']
yNDpfin = inf['NDp']
yRHfin = inf['RH']
# yRHefin = inf['RH_e']
# yRHhfin = inf['RH_h']
# yFHefin = inf['FHall_e']
# yFHhfin = inf['FHall_h']
# print("")
if cparams.print_level >= 1:
print("{:8}\t{:8}\t{:10}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}"
.format('T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,fin', 'Nh,fin', 'NA-,fin', 'ND+,fin', 'RH,fin', '|Ns,fin|'))
for i in range(len(dos.xT)):
print("{:8.4f}\t{:8.4f}\t{:10.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}"
.format(dos.xT[i], xTinv[i], yEFfin[i], dos.yN[i], yNefin[i], yNhfin[i], yNAmfin[i], yNDpfin[i], yRHfin[i], yNsabsfin[i]))
if cparams.mode == 'fit_RH':
yNc = yRHfin
if cparams.mode == 'fit_h':
yNc = yNhfin
elif cparams.mode == 'fit_e':
yNc = yNefin
else:
yNc = yNsabsfin
S2 = 0.0
for iT in range(len(dos.xT)):
T = dos.xT[iT]
if cparams.mode == 'fit_RH':
if dos.yRH[iT] < 0.0:
lnRH = -log(-dos.yRH[iT])
else:
lnRH = log(dos.yRH[iT]+eps)
if yRHfin[iT] < 0.0:
lnRHfin = -log(-yRHfin[iT])
else:
lnRHfin = log(yRHfin[iT]+eps)
# S2 += (log(dos.yRH[iT]+eps) - log(yRHfin[iT]+eps))**2
S2 += (lnRH - lnRHfin)**2
else:
S2 += (log(dos.yN[iT]+eps) - log(yNc[iT]+eps))**2
S2 /= (len(dos.xT) - 1)
dos.set_parameters(aiorg)
return S2
[ドキュメント]
def diff1(app, cparams, dos, mu, ai):
"""
概要: 二乗誤差関数 `calS2` の一次導関数を数値的に近似計算します。
詳細説明:
各最適化パラメータ `ai` に対して微小な変化 `h_diff` を与え、
前方差分法 (forward difference method) を用いて `calS2` の勾配を近似します。
この導関数は、共役勾配法 (CG) など一部の最適化アルゴリズムで利用されます。
:param app: アプリケーションオブジェクト。
:type app: tkApplication
:param cparams: アプリケーションパラメータオブジェクト。
:type cparams: tkParams
:param dos: 状態密度オブジェクト。
:type dos: tkDOS
:param mu: 移動度オブジェクト。
:type mu: tkMobility
:param ai: 現在の最適化パラメータ値のリスト。
:type ai: list[float]
:returns: 各パラメータに関する一次導関数の配列。
:rtype: numpy.ndarray
"""
n = len(ai)
f0 = calS2(app, cparams, dos, mu, ai) # Pass app, cparams, dos, mu to calS2
df = np.empty(n)
for i in range(n):
aii = ai.copy()
aii[i] = ai[i] + cparams.h_diff # Use cparams.h_diff
df[i] = (calS2(app, cparams, dos, mu, aii) - f0) / cparams.h_diff # Pass app, cparams, dos, mu to calS2
return df
[ドキュメント]
def plot(app, cparams, dos, mobility, xTinv, yN, yNhini, yNeini, yNsabsini, yNhfin = None, yNefin = None, yNsabsfin = None,
ymu = None, ymu_fin = None, yNhHallfin = None, yNeHallfin = None):
"""
概要: キャリア濃度、移動度、ホール係数などの温度依存性データをプロットします。
詳細説明:
観測データと初期計算値、最終フィッティング値を比較してグラフを作成します。
通常は2x2のサブプロット構成で、キャリア濃度（対数スケール）、
フィッティング後のキャリア濃度、移動度を表示します。
必要に応じてホールキャリア濃度などの詳細なフィッティング結果も表示します。
:param app: アプリケーションオブジェクト。
:type app: tkApplication
:param cparams: アプリケーションパラメータオブジェクト。
:type cparams: tkParams
:param dos: 状態密度オブジェクト。
:type dos: tkDOS
:param mobility: 移動度オブジェクト。
:type mobility: tkMobility
:param xTinv: 温度の逆数 (1000/T) のリスト。
:type xTinv: list[float]
:param yN: 観測されたキャリア濃度 (cm^-3) のリスト。
:type yN: list[float]
:param yNhini: 初期計算されたホール濃度 (cm^-3) のリスト。
:type yNhini: list[float]
:param yNeini: 初期計算された電子濃度 (cm^-3) のリスト。
:type yNeini: list[float]
:param yNsabsini: 初期計算された絶対キャリア濃度 (|Ns|) (cm^-3) のリスト。
:type yNsabsini: list[float]
:param yNhfin: (オプション) 最終フィッティングされたホール濃度 (cm^-3) のリスト。
:type yNhfin: list[float] or None
:param yNefin: (オプション) 最終フィッティングされた電子濃度 (cm^-3) のリスト。
:type yNefin: list[float] or None
:param yNsabsfin: (オプション) 最終フィッティングされた絶対キャリア濃度 (|Ns|) (cm^-3) のリスト。
:type yNsabsfin: list[float] or None
:param ymu: (オプション) 観測された移動度 (cm^2/(Vs)) のリスト。
:type ymu: list[float] or None
:param ymu_fin: (オプション) 最終フィッティングされた移動度 (cm^2/(Vs)) のリスト。
:type ymu_fin: list[float] or None
:param yNhHallfin: (オプション) 最終フィッティングされたホールによるホール濃度 (cm^-3) のリスト。
:type yNhHallfin: list[float] or None
:param yNeHallfin: (オプション) 最終フィッティングされた電子によるホール濃度 (cm^-3) のリスト。
:type yNeHallfin: list[float] or None
:returns: None
:rtype: None
"""
kplot = [0.01, 2.0]
plt = app.plt
fig = app.fig
fontsize = cparams.fontsize
legend_fontsize = cparams.legend_fontsize
plt.clf()
ax1 = fig.add_subplot(2, 2, 1)
ax2 = fig.add_subplot(2, 2, 2)
axmu = fig.add_subplot(2, 2, 3)
plt.title(cparams.infile, fontsize = fontsize)
ax1.plot(xTinv, yN, label = 'N(obs)', linestyle = 'none', marker = 'o')
ax1.plot(xTinv, yNhini, label = '$N_h$(ini)', linestyle = 'dashed', color = 'red', linewidth = 0.5)
ax1.plot(xTinv, yNeini, label = '$N_e$(ini)', linestyle = 'dashed', color = 'blue', linewidth = 0.5)
ax1.plot(xTinv, yNsabsini, label = '|$N_s$(ini)|', linestyle = '-', color = 'green', linewidth = 0.5)
ax1.set_yscale('log')
ax1.set_ylim([min(yN) * kplot[0], max(yN) * kplot[1]])
ax1.set_xlabel('1000/T (K$^{-1}$)', fontsize = fontsize)
ax1.set_ylabel('N (cm$^{-3}$)', fontsize = fontsize)
ax1.legend(fontsize = legend_fontsize)
ax1.tick_params(labelsize = fontsize)
if yNhfin:
print_parameters(dos.varname, dos.varunit, dos.parameter_list(), dos.optid)
ax2.plot(xTinv, yN, label = '$N$(obs,Hall)', linestyle = 'none', marker = 'o')
ax2.plot(xTinv, yNhfin, label = '$N_h$(fin)', linestyle = 'dashed', color = 'red', linewidth = 0.5)
ax2.plot(xTinv, yNefin, label = '$N_e$(fin)', linestyle = 'dashed', color = 'blue', linewidth = 0.5)
ax2.plot(xTinv, yNsabsfin, label = '|$N_s$(fin)|', linestyle = '-', color = 'green', linewidth = 0.5)
if yNhHallfin:
ax2.plot(xTinv, yNhHallfin, label = '$N_{h,Hall}$(fin)', linestyle = 'dotted', color = 'red', linewidth = 0.3)
ax2.plot(xTinv, yNeHallfin, label = '$N_{e,Hall}$(fin)', linestyle = 'dotted', color = 'blue', linewidth = 0.3)
if yNhfin:
ax2.set_yscale('log')
ax2.set_ylim([min(yN) * kplot[0], max(yN) * kplot[1]])
ax2.set_xlabel('1000/T (K$^{-1}$)', fontsize = fontsize)
ax2.set_ylabel('N (cm$^{-3}$)', fontsize = fontsize)
ax2.legend(fontsize = legend_fontsize)
ax2.tick_params(labelsize = fontsize)
if ymu is not None:
axmu.plot(dos.xT, dos.ymu, label = '$\mu$(obs,Hall)', color = 'black', linewidth = 0.5, linestyle = '-', marker = 'o')
if ymu_fin is not None:
axmu.plot(dos.xT, ymu_fin, label = '$\mu$(fin)', color = 'blue', linewidth = 0.5, linestyle = 'dashed', marker = '^')
axmu.set_xlabel("T (K)", fontsize = fontsize)
axmu.set_ylabel("$\mu$ (cm$^2/(Vs)$)", fontsize = fontsize)
axmu.legend(fontsize = legend_fontsize)
axmu.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.01)
[ドキュメント]
def callback(app, cparams, dos, mobility, xk):
"""
概要: `scipy.optimize.minimize` のコールバック関数として、最適化の進捗状況を更新・表示します。
詳細説明:
最適化の各イテレーションで呼び出され、現在のパラメータセット `xk` を用いて
二乗誤差 `S2` を計算します。`outputinterval` で指定された間隔で、
現在のパラメータと `S2` の値を標準出力に表示します。
また、`graphupdateinterval` で指定された間隔でグラフを更新し、
フィッティングの進行状況を視覚的に確認できるようにします。
:param app: アプリケーションオブジェクト。
:type app: tkApplication
:param cparams: アプリケーションパラメータオブジェクト。
:type cparams: tkParams
:param dos: 状態密度オブジェクト。
:type dos: tkDOS
:param mobility: 移動度オブジェクト。
:type mobility: tkMobility
:param xk: 現在の最適化パラメータの配列。
:type xk: numpy.ndarray
:returns: None
:rtype: None
"""
S2 = calS2(app, cparams, dos, mobility, xk)
ai0 = dos.parameter_list()
ainew = recover_parameters(xk, dos.optid, ai0)
dos.set_parameters(ainew)
if app.iter % cparams.outputinterval == 0:
print("callback {:04d}: ".format(app.iter), end = '')
for i in range(len(ainew)):
print("{:10.4g} ".format(ainew[i]), end = '')
print(" S2={:10.6g}".format(S2))
if app.iter % cparams.graphupdateinterval == 0:
ymufin = []
xTinv, yEFfin, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, yNsabsfin, ydQ \
= dos.cal_T_lists(dos.xT, print_level = cparams.print_level)
ymu_fit = []
for i in range(len(dos.xT)):
ymu_fit.append(dos.ys[i] / e / (yNefin[i] + yNhfin[i]))
plot(app, cparams, dos, mobility, xTinv, dos.yN, dos.yNhini, dos.yNeini, dos.yNsabsini,
yNhfin, yNefin, yNsabsfin, ymu = dos.ymu, ymu_fin = ymu_fit)
app.iter += 1
try:
pass
except KeyboardInterrupt:
print ('KeyboardInterrupt exception is caught')
terminate(usage = usage)
[ドキュメント]
def callbackRH(app, cparams, dos, mobility, xk):
"""
概要: ホール係数フィッティングモードにおける `scipy.optimize.minimize` のコールバック関数です。
詳細説明:
`callback` 関数と同様の機能ですが、`calS2RH` を使用して二乗誤差を計算し、
`fit_RH` モードに特化した情報（ホールキャリア濃度など）を含むグラフを更新します。
最適化の進行状況とホール係数のフィッティング状況を可視化します。
:param app: アプリケーションオブジェクト。
:type app: tkApplication
:param cparams: アプリケーションパラメータオブジェクト。
:type cparams: tkParams
:param dos: 状態密度オブジェクト。
:type dos: tkDOS
:param mobility: 移動度オブジェクト。
:type mobility: tkMobility
:param xk: 現在の最適化パラメータの配列。
:type xk: numpy.ndarray
:returns: None
:rtype: None
"""
S2 = calS2RH(app, cparams, dos, mobility, xk)
ai0 = dos.parameter_list()
ainew = recover_parameters(xk, dos.optid, ai0)
dos.set_parameters(ainew)
if app.iter % cparams.outputinterval == 0:
print("callback {:04d}: ".format(app.iter), end = '')
for i in range(len(ainew)):
print("{:10.4g} ".format(ainew[i]), end = '')
print(" S2={:10.6g}".format(S2))
if app.iter % cparams.graphupdateinterval == 0:
ymufin = []
# xTinv, yEFfin, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, yNsabsfin, ydQ = dos.cal_T_lists(dos.xT, print_level = cparams.print_level)
inf = dos.cal_T_lists2(dos.xT, mobility, print_level = cparams.print_level)
xTinv = inf['Tinv']
# yEFfin = inf['EF']
yNefin = inf['n_e']
yNhfin = inf['n_h']
yNeHallfin = inf['nHall_e']
yNhHallfin = inf['nHall_h']
yNsfin = inf['n_s']
yNsabsfin = inf['n_s(abs)']
# yNAmfin = inf['NAm']
# yNDpfin = inf['NDp']
# yRHfin = inf['RH']
# yRHefin = inf['RH_e']
# yRHhfin = inf['RH_h']
# yFHefin = inf['FHall_e']
# yFHhfin = inf['FHall_h']
ymu_fit = []
for i in range(len(dos.xT)):
ymu_fit.append(dos.ys[i] / e / (yNefin[i] + yNhfin[i]))
plot(app, cparams, dos, mobility, dos.xTinv, dos.yN, dos.yNhini, dos.yNeini, dos.yNsabsini,
yNhfin, yNefin, yNsabsfin, ymu = dos.ymu, ymu_fin = ymu_fit,
yNhHallfin = yNhHallfin, yNeHallfin = yNeHallfin)
app.iter += 1
try:
pass
except KeyboardInterrupt:
print ('KeyboardInterrupt exception is caught')
terminate(usage = usage)
[ドキュメント]
def fit(app, cparams, dos, mobility):
"""
概要: 温度依存性ホールキャリア濃度データに対して非線形最小二乗フィッティングを行います。
詳細説明:
ホール効果測定によって得られたキャリア濃度データ (Ns = 1/eRH) を読み込み、
`cparams.mode` に応じて単一キャリアまたは二キャリアモデル（散乱因子 `r=0.5` 固定）
を用いてDOSパラメータの最適化を行います。
最適化には `scipy.optimize.minimize` が使用され、指定されたアルゴリズム (`cparams.method`)
に基づいてパラメータが調整されます。
最適化されたパラメータはファイルに保存され、初期値とフィッティング後の結果がプロットされます。
:param app: アプリケーションオブジェクト。
:type app: tkApplication
:param cparams: アプリケーションパラメータオブジェクト。
:type cparams: tkParams
:param dos: 状態密度オブジェクト。
:type dos: tkDOS
:param mobility: 移動度オブジェクト。
:type mobility: tkMobility
:returns: None
:rtype: None
"""
print("Fitting to mu-T Hall data by nonlinear least-squares method")
print(" r: {}".format(mobility.rfac))
# read data
print("")
print("Read Hall data from {}".format(cparams.infile))
header, xT, yN, ymu, ys = dos.read_excel_Hall(cparams.infile, usage = usage)
# header, xT, yN, ymu, ys, yRH = dos.read_excel_Hall2(cparams.infile, usage = usage)
print(" Headers:", header)
if xT is None:
terminate(f"Temperature data is not found in [{cparams.infile}]", usage = usage)
if yN is None:
terminate(f"Carrier density data is not found in [{cparams.infile}]", usage = usage)
if ymu is None:
terminate(f"Mobility data is not found in [{cparams.infile}]", usage = usage)
if ys is None:
terminate(f"Conductivity data is not found in [{cparams.infile}]", usage = usage)
# calculate initial values
xTinv = []
ymuini = []
print("")
print("Calculate mu-T from the initial parameters")
# print("{:8}\t{:12}".format("T(K)", "mu(cm2/Vs)"))
dos.xTinv, dos.yEFini, dos.yNeini, dos.yNhini, dos.yNAmini, dos.yNDpini, dos.yRHini, dos.yNsini, dos.yNsabsini, dos.ydQ \
= dos.cal_T_lists(xT, print_level = cparams.print_level)
print("")
print("{:8}\t{:8}\t{:10}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}"
.format('T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,fin', 'Nh,fin', 'NA-,fin', 'ND+,fin', 'RH,fin', 'Ns,fin'))
for i in range(len(xT)):
print("{:8.4f}\t{:8.4f}\t{:10.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}"
.format(dos.xT[i], dos.xTinv[i], dos.yEFini[i], dos.yN[i],
dos.yNeini[i], dos.yNhini[i], dos.yNAmini[i], dos.yNDpini[i], dos.yRHini[i], dos.yNsini[i]))
# exit()
#=============================
# グラフの表示
#=============================
print("")
print("plot")
app.plt = plt
app.fig = plt.figure(figsize = cparams.figsize)
plot(app, cparams, dos, mobility, dos.xTinv, dos.yN, dos.yNhini, dos.yNeini, dos.yNsabsini)
#=============================
# Optimization
#=============================
print("")
print("")
print("Nonlinear least-squares fitting by ", cparams.method)
print(" tol=", cparams.tol)
app.iter = 0
ai0 = dos.parameter_list()
ai = choose_parameters(ai0, dos.optid)
ret = minimize(lambda xi: calS2(app, cparams, dos, mobility, xi), ai, method = cparams.method,
jac = lambda xi: diff1(app, cparams, dos, mobility, xi),
tol = cparams.tol, callback = lambda xk: callback(app, cparams, dos, mobility, xk),
options = {'maxiter': cparams.maxiter, "disp": True})
# print("ret=", ret)
if cparams.method == 'nelder-mead':
simplex = ret['final_simplex']
ai = simplex[0][0]
fmin = ret['fun']
elif cparams.method == 'cg':
ai = ret['x']
fmin = ret['fun']
elif cparams.method == 'powell':
ai = ret['x']
fmin = ret['fun']
elif cparams.method == 'bfgs':
ai = ret['x']
fmin = calS2(app, cparams, dos, mobility, ai) # Use calS2 with correct args
aifin = recover_parameters(ai, dos.optid, ai0)
dos.set_parameters(aifin)
print("")
print("Optimized at S2={:12.6g}:".format(fmin))
print_parameters(dos.varname, dos.varunit, dos.parameter_list(), dos.optid)
dos.NC = dos.meff2NC_FEA(dos.meeff, cparams.T0)
dos.NV = dos.meff2NC_FEA(dos.mheff, cparams.T0)
dos.DC0 = dos.meff2DC0_FEA(dos.meeff, cparams.T0)
dos.DV0 = dos.meff2DC0_FEA(dos.mheff, cparams.T0)
print(" NV : {} me at {} K".format(dos.NV, cparams.T0))
print(" DV0: {} me at {} K".format(dos.DV0, cparams.T0))
print(" NC : {} me at {} K".format(dos.NC, cparams.T0))
print(" DC0: {} me at {} K".format(dos.DC0, cparams.T0))
print("")
print("Save to [{}]".format(cparams.parameterfile))
cparams.save_parameters(cparams.parameterfile, "Preferences")
dos.save_parameters(cparams.parameterfile, dos.optid,
{"last_fitmode": cparams.mode, "rfac": mobility.rfac, "l0": mobility.l0, 'S2N': fmin})
print("")
print("Calculate final data")
xTinv, yEFfin, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, yNsabsfin, ydQ \
= dos.cal_T_lists(xT, print_level = cparams.print_level)
print("")
print("{:8}\t{:8}\t{:10}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}"
.format('T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,fin', 'Nh,fin', 'NA-,fin', 'ND+,fin', 'RH,fin', 'Ns,fin', 'mu(Hall)', 'mu,fin'))
ymu_fit = []
for i in range(len(xT)):
ymu_fit.append(ys[i] / e / (yNefin[i] + yNhfin[i]))
print("{:8.4f}\t{:8.4f}\t{:10.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}"
.format(xT[i], xTinv[i], yEFfin[i], yN[i], yNefin[i], yNhfin[i], yNAmfin[i], yNDpfin[i], yRHfin[i], yNsfin[i], ymu[i], ymu_fit[i]))
print("Save fitting data to [{}]".format(cparams.outcsvfile))
save_csv(cparams.outcsvfile,
['T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,fin', 'Nh,fin', 'NA-,fin', 'ND+,fin', 'RH,fin', 'Ns,fin', 'mu(Hall)(cm2/Vs)', 'mu,fin'],
[xT, xTinv, yEFfin, yN, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, ymu, ymu_fit])
#=============================
# グラフの表示
#=============================
print("")
print("Plot optimized")
plot(app, cparams, dos, mobility, dos.xTinv, dos.yN, dos.yNhini, dos.yNeini, dos.yNsabsini,
yNhfin, yNefin, yNsabsfin, ymu = ymu, ymu_fin = ymu_fit)
print("")
print("Press ENTER to exit>>", end = '')
input()
print("")
app.terminate(usage = usage)
[ドキュメント]
def fit_RH(app, cparams, dos, mobility):
"""
概要: ホール係数 (RH) の温度依存性データに対して非線形最小二乗フィッティングを行います。
詳細説明:
ホール効果測定によって得られたホール係数データ (`yRH`) を読み込み、
二キャリアモデルと、`mobility.rfac` で設定された散乱因子を用いて、
DOSパラメータの最適化を行います。
最適化には `scipy.optimize.minimize` が使用され、指定されたアルゴリズム
(`cparams.method`) に基づいてパラメータが調整されます。
最適化されたパラメータはファイルに保存され、初期値とフィッティング後の結果が
詳細な情報（ホールホール濃度、電子ホール濃度など）とともにプロットされます。
:param app: アプリケーションオブジェクト。
:type app: tkApplication
:param cparams: アプリケーションパラメータオブジェクト。
:type cparams: tkParams
:param dos: 状態密度オブジェクト。
:type dos: tkDOS
:param mobility: 移動度オブジェクト。
:type mobility: tkMobility
:returns: None
:rtype: None
"""
print("Fitting to mu-T Hall data by nonlinear least-squares method")
print(" r: {}".format(mobility.rfac))
print(" Calculate two-carrier model RH:")
print(f" mu_h: {dos.cal_mobility_h_description}")
print(f" mu_e: {dos.cal_mobility_e_description}")
# read data
print("")
print("Read Hall data from {}".format(cparams.infile))
header, xT, yN, ymu, ys, yRH = dos.read_excel_Hall2(cparams.infile, usage = usage)
print(" Headers:", header)
# print(" xT :", xT)
# print(" yN :", yN)
# print(" ymu:", ymu)
# print(" ys :", ys)
# print(" yRH:", yRH)
if xT is None:
app.terminate(f"Temperature data is not found in [{cparams.infile}]", usage = usage)
if yN is None:
app.terminate(f"Carrier density data is not found in [{cparams.infile}]", usage = usage)
if ymu is None:
app.terminate(f"Mobility data is not found in [{cparams.infile}]", usage = usage)
if ys is None:
app.terminate(f"Conductivity data is not found in [{cparams.infile}]", usage = usage)
if yRH is None:
app.terminate(f"RH data is not found in [{cparams.infile}]", usage = usage)
# calculate initial values
xTinv = []
ymuini = []
print("")
print("Calculate mu-T from the initial parameters")
# print("{:8}\t{:12}".format("T(K)", "mu(cm2/Vs)"))
# xTinv, yEFini, yNeini, yNhini, yNAmini, yNDpini, yRHini, yNsini, yNsabsini, ydQ = dos.cal_T_lists(xT, print_level = cparams.print_level)
inf = dos.cal_T_lists2(xT, mobility, print_level = cparams.print_level)
dos.xTinv = inf['Tinv']
dos.yEFini = inf['EF']
dos.yNeini = inf['n_e']
dos.yNhini = inf['n_h']
dos.yNeHallini = inf['nHall_e']
dos.yNhHallini = inf['nHall_h']
dos.yNsini = inf['n_s']
dos.yNsabsini = inf['n_s(abs)']
dos.yNAmini = inf['NAm']
dos.yNDpini = inf['NDp']
dos.yRHini = inf['RH']
dos.yRHeini = inf['RH_e']
dos.yRHhini = inf['RH_h']
dos.yFHeini = inf['FHall_e']
dos.yFHhini = inf['FHall_h']
print("")
print("{:8}\t{:8}\t{:10}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}"
.format('T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,ini', 'Nh,ini', 'NA-,ini', 'ND+,ini', 'RH,ini', 'Ns,ini'))
for i in range(len(xT)):
print("{:8.4f}\t{:8.4f}\t{:10.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}"
.format(dos.xT[i], dos.xTinv[i], dos.yEFini[i], dos.yN[i],
dos.yNeini[i], dos.yNhini[i], dos.yNAmini[i], dos.yNDpini[i], dos.yRHini[i], dos.yNsini[i]))
# exit()
#=============================
# グラフの表示
#=============================
print("")
print("plot")
app.plt = plt
app.fig = plt.figure(figsize = cparams.figsize)
plot(app, cparams, dos, mobility, dos.xTinv, dos.yN, dos.yNhini, dos.yNeini, dos.yNsabsini)
#=============================
# Optimization
#=============================
print("")
print("")
print("Nonlinear least-squares fitting by ", cparams.method)
print(" tol=", cparams.tol)
app.iter = 0
ai0 = dos.parameter_list()
ai = choose_parameters(ai0, dos.optid)
ret = minimize(lambda xk: calS2RH(app, cparams, dos, mobility, xk), ai, method = cparams.method,
jac = lambda xk: diff1(app, cparams, dos, mobility, xk),
tol = cparams.tol,
callback = lambda xk: callbackRH(app, cparams, dos, mobility, xk),
options = {'maxiter': cparams.maxiter, "disp":True})
# print("ret=", ret)
if cparams.method == 'nelder-mead':
simplex = ret['final_simplex']
ai = simplex[0][0]
fmin = ret['fun']
elif cparams.method == 'cg':
ai = ret['x']
fmin = ret['fun']
elif cparams.method == 'powell':
ai = ret['x']
fmin = ret['fun']
elif cparams.method == 'bfgs':
ai = ret['x']
fmin = calS2RH(app, cparams, dos, mobility, ai) # Use calS2RH with correct args
aifin = recover_parameters(ai, dos.optid, ai0)
dos.set_parameters(aifin)
print("")
print("Optimized at S2={:12.6g}:".format(fmin))
print_parameters(dos.varname, dos.varunit, dos.parameter_list(), dos.optid)
dos.NC = dos.meff2NC_FEA(dos.meeff, cparams.T0)
dos.NV = dos.meff2NC_FEA(dos.mheff, cparams.T0)
dos.DC0 = dos.meff2DC0_FEA(dos.meeff, cparams.T0)
dos.DV0 = dos.meff2DC0_FEA(dos.mheff, cparams.T0)
print(" NV : {} me at {} K".format(dos.NV, cparams.T0))
print(" DV0: {} me at {} K".format(dos.DV0, cparams.T0))
print(" NC : {} me at {} K".format(dos.NC, cparams.T0))
print(" DC0: {} me at {} K".format(dos.DC0, cparams.T0))
print("")
print("Save to [{}]".format(cparams.parameterfile))
cparams.save_parameters(cparams.parameterfile, "Preferences")
dos.save_parameters(cparams.parameterfile, dos.optid,
{"last_fitmode": cparams.mode, "rfac": mobility.rfac, "l0": mobility.l0, 'S2N': fmin})
print("")
print("Calculate final data")
# xTinv, yEFfin, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, yNsabsfin, ydQ = dos.cal_T_lists(xT, print_level = cparams.print_level)
inf = dos.cal_T_lists2(xT, mobility, print_level = cparams.print_level)
xTinv = inf['Tinv']
yEFfin = inf['EF']
yNefin = inf['n_e']
yNhfin = inf['n_h']
yNeHallfin = inf['nHall_e']
yNhHallfin = inf['nHall_h']
yNsfin = inf['n_s']
yNsabsfin = inf['n_s(abs)']
yNAmfin = inf['NAm']
yNDpfin = inf['NDp']
yRHfin = inf['RH']
yRHefin = inf['RH_e']
yRHhfin = inf['RH_h']
yFHefin = inf['FHall_e']
yFHhfin = inf['FHall_h']
print("")
print("{:8}\t{:8}\t{:10}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}"
.format('T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,fin', 'Nh,fin', 'NA-,fin', 'ND+,fin', 'RH,fin', 'Ns,fin', 'mu(Hall)', 'mu,fin'))
ymu_fit = []
for i in range(len(xT)):
ymu_fit.append(ys[i] / e / (yNefin[i] + yNhfin[i]))
print("{:8.4f}\t{:8.4f}\t{:10.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}"
.format(xT[i], xTinv[i], yEFfin[i], yN[i], yNefin[i], yNhfin[i], yNAmfin[i], yNDpfin[i], yRHfin[i], yNsfin[i], ymu[i], ymu_fit[i]))
print("Save fitting data to [{}]".format(cparams.outcsvfile))
save_csv(cparams.outcsvfile,
['T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,fin', 'Nh,fin', 'NA-,fin', 'ND+,fin', 'RH,fin', 'Ns,fin', 'mu(Hall)(cm2/Vs)', 'mu,fin',
'FH,e', 'FH,h', 'Ne,Hall', 'Nh,Hall'],
[xT, xTinv, yEFfin, yN, yNefin, yNhfin, yNAmfin, yNDpfin, yRHfin, yNsfin, ymu, ymu_fit,
yFHefin, yFHhfin, yNeHallfin, yNhHallfin])
#=============================
# グラフの表示
#=============================
print("")
print("Plot optimized")
plot(app, cparams, dos, mobility, dos.xTinv, dos.yN, dos.yNhini, dos.yNeini, dos.yNsabsini,
yNhfin, yNefin, yNsabsfin, ymu = dos.ymu, ymu_fin = ymu_fit,
yNhHallfin = yNhHallfin, yNeHallfin = yNeHallfin)
print("")
print("Press ENTER to exit>>", end = '')
input()
print("")
app.terminate(usage = usage)
[ドキュメント]
def sim(app, cparams, dos, mobility):
"""
概要: 現在のパラメータに基づいてホール効果関連データの温度依存性をシミュレーションし、プロットします。
詳細説明:
入力ファイルから観測データを読み込み、現在のDOSおよび移動度パラメータを用いて、
キャリア濃度、移動度、ホール係数、活性化エネルギーなどの温度依存性を計算します。
計算結果は標準出力に表示され、CSVファイルに保存されます。
また、キャリア濃度、移動度、ホール因子、フェルミ準位、活性化エネルギーなどを含む
複数のグラフが生成され、視覚的に結果を確認できます。
:param app: アプリケーションオブジェクト。
:type app: tkApplication
:param cparams: アプリケーションパラメータオブジェクト。
:type cparams: tkParams
:param dos: 状態密度オブジェクト。
:type dos: tkDOS
:param mobility: 移動度オブジェクト。
:type mobility: tkMobility
:returns: None
:rtype: None
"""
print("")
print("Simulate T dependences")
print(" r: {}".format(mobility.rfac))
# read data
print("")
print("Read Hall data from {}".format(cparams.infile))
# header, xT, yN, ymu, ys = dos.read_excel_Hall(cparams.infile, usage = usage)
header, xT, yN, ymu, ys, yRH = dos.read_excel_Hall2(cparams.infile, usage = usage)
print(" Headers:", header)
print(" xT :", xT)
print(" yN :", yN)
print(" ymu:", ymu)
print(" ys :", ys)
print(" yRH:", yRH)
# calculate initial values
print("")
# xTinv, yEFini, yNeini, yNhini, yNAmini, yNDpini, yRHini, yNsini, yNsabsini, ydQ = dos.cal_T_lists(xT, print_level = cparams.print_level)
inf = dos.cal_T_lists2(xT, mobility, print_level = cparams.print_level)
xTinv = inf['Tinv']
yEFini = inf['EF']
yNeini = inf['n_e']
yNhini = inf['n_h']
yNeHallini = inf['nHall_e']
yNhHallini = inf['nHall_h']
yNsini = inf['n_s']
yNsabsini = inf['n_s(abs)']
yNAmini = inf['NAm']
yNDpini = inf['NDp']
yRHini = inf['RH']
yRHeini = inf['RH_e']
yRHhini = inf['RH_h']
yFHeini = inf['FHall_e']
yFHhini = inf['FHall_h']
yNT23 = []
yNT43 = []
print("Calculate n-T from the initial parameters")
print("{:8}\t{:8}\t{:10}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}\t{:14}"
.format('T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,ini', 'Nh,ini', 'NA-,ini', 'ND+,ini', 'RH,ini', 'Ns,ini', 'mu(Hall)', 'mu,ini'))
ymu_ini = []
for i in range(len(xT)):
print("i=", i, ys[i], yNeini[i], yNhini[i])
ymu_ini.append(ys[i] / e / (yNeini[i] + yNhini[i]))
yNT23.append(yN[i] * pow(xT[i], -2.0 / 3.0))
yNT43.append(yN[i] * pow(xT[i], -4.0 / 3.0))
print("{:8.4f}\t{:8.4f}\t{:10.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}\t{:14.6g}"
.format(xT[i], xTinv[i], yEFini[i], yN[i], yNeini[i], yNhini[i], yNAmini[i], yNDpini[i],
yRHini[i], yNsini[i], ymu[i], ymu_ini[i]))
print("")
# xTinv, yEFini, yNeini, yNhini, yNAmini, yNDpini, yRHini, yNsini, yNsabsini, ydQ = dos.cal_T_lists(xT, print_level = dos.print_level)
Tinv, yEa = Ea_Arrhenius(xT, yN, eps = 1.0e-300, kTinv = 1000.0)
Tinv, yEaT23 = Ea_Arrhenius(xT, yNT23, eps = 1.0e-300, kTinv = 1000.0)
Tinv, yEaT43 = Ea_Arrhenius(xT, yNT43, eps = 1.0e-300, kTinv = 1000.0)
nT = len(xT)
print("Activation energy")
print("{:8}\t{:14}\t{:14}\t{:14}".format("T(K)", "Ea(eV)(log(N))", "Ea(eV)(log(N*T^-2/3))", "Ea(eV)(log(N*T^-4/3))"))
for iT in range(nT):
print("{:8.4g}\t{:14.6g}\t{:14.6g}\t{:14.6g}".format(xT[iT], yEa[iT], yEaT23[iT], yEaT43[iT]))
print("")
print("Save to [{}]".format(cparams.parameterfile))
cparams.save_parameters(cparams.parameterfile, "Preferences")
dos.save_parameters(cparams.parameterfile, dos.optid, {"rfac": mobility.rfac, "l0": mobility.l0})
print("Save simulation data to [{}]".format(cparams.outcsvfile))
save_csv(cparams.outcsvfile,
['T', '1000/T', 'EF(eV)', 'Ne,obs(cm-3)', 'Ne,ini', 'Nh,ini', 'NA-,ini', 'ND+,ini', 'RH,ini', 'Ns,ini', 'mu(Hall)(cm2/Vs)', 'mu,ini',
'FH,e', 'FH,h', 'Ne,Hall', 'Nh,Hall'],
[xT, xTinv, yEFini, yN, yNeini, yNhini, yNAmini, yNDpini, yRHini, yNsini, ymu, ymu_ini,
yFHeini, yFHhini, yNeHallini, yNhHallini])
#=============================
# グラフの表示
#=============================
app.fig = plt.figure(figsize = cparams.figsize)
fig = app.fig
fontsize = cparams.fontsize
legend_fontsize = cparams.legend_fontsize
ax1 = fig.add_subplot(2, 3, 1)
ax3 = fig.add_subplot(2, 3, 2)
axmu = fig.add_subplot(2, 3, 3)
ax2 = fig.add_subplot(2, 3, 4)
ax4 = fig.add_subplot(2, 3, 5)
axFH = fig.add_subplot(2, 3, 6)
axEF = axFH.twinx()
ax1.set_title(cparams.infile, fontsize = fontsize)
ins1 = axFH.plot(xTinv, yFHeini, label = '$F_{Hall,e}$', color = 'blue', linestyle = '-', linewidth = 0.5)
ins2 = axFH.plot(xTinv, yFHhini, label = '$F_{Hall,h}$', color = 'red', linestyle = '-', linewidth = 0.5)
ins3 = axEF.plot(xTinv, yEFini, label = '$E_F$', color = 'black', linestyle = '-', linewidth = 0.5)
ins4 = axEF.plot(axFH.get_xlim(), [dos.EV, dos.EV], label = '$E_V$', color = 'green', linestyle = 'dashed', linewidth = 0.3)
ins5 = axEF.plot(axFH.get_xlim(), [dos.EC, dos.EC], label = '$E_C$', color = 'green', linestyle = 'dashed', linewidth = 0.3)
axFH.set_xlabel("1000/T (K$^{-1}$)", fontsize = fontsize)
axFH.set_ylabel("$F_{Hall}$ (cm$^{-3}$/C)", fontsize = fontsize)
axEF.set_ylabel("$E_F$ (eV)", fontsize = fontsize)
ins = ins1 + ins2 + ins3 + ins4 + ins5
axFH.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
# axFH.legend(fontsize = legend_fontsize)
axFH.tick_params(labelsize = fontsize)
ax1.plot(xTinv, yN, label = '$N(obs)$', color = 'black', linestyle = '', marker = 'o')
ax1.plot(xTinv, yNeini, label = '$N_e(ini)$', color = 'blue', linewidth = 0.5, linestyle = '-')
ax1.plot(xTinv, yNhini, label = '$N_h(ini)$', color = 'red', linewidth = 0.5, linestyle = '-')
ax1.plot(xTinv, yNeHallini, label = '$N_{e,Hall}(ini)$', color = 'blue', linewidth = 0.5, linestyle = 'dashed')
ax1.plot(xTinv, yNhHallini, label = '$N_{h,Hall}(ini)$', color = 'red', linewidth = 0.5, linestyle = 'dashed')
ax1.plot(xTinv, yNsabsini, label = '$N_s(ini)$', color = 'purple', linewidth = 0.5, linestyle = 'dashed')
ax1.set_xlabel("1000/T (K$^{-1}$)", fontsize = fontsize)
ax1.set_ylabel("N (cm$^{-3}$)", fontsize = fontsize)
ax1.set_yscale('log')
ax1.legend(fontsize = legend_fontsize)
ax1.tick_params(labelsize = fontsize)
ax3.plot(xT, yN, label = '$N(obs)$', color = 'black', linestyle = '', marker = 'o')
ax3.plot(xT, yNeini, label = '$N_e(ini)$', color = 'blue', linewidth = 0.5, linestyle = '-')
ax3.plot(xT, yNhini, label = '$N_h(ini)$', color = 'red', linewidth = 0.5, linestyle = '-')
ax3.plot(xT, yNsabsini, label = '$N_s(ini)$', color = 'purple', linewidth = 0.5, linestyle = 'dashed')
ax3.set_xlabel("T (K)", fontsize = fontsize)
ax3.set_ylabel("N (cm$^{-3}$)", fontsize = fontsize)
ax3.set_yscale('log')
ax3.legend(fontsize = legend_fontsize)
ax3.tick_params(labelsize = fontsize)
axmu.plot(xT, ymu, label = '$\mu$(Hall)', color = 'black', linewidth = 0.5, linestyle = '-', marker = 'o')
axmu.plot(xT, ymu_ini, label = '$\mu$(ini)', color = 'blue', linewidth = 0.5, linestyle = 'dashed', marker = '^')
axmu.set_xlabel("T (K)", fontsize = fontsize)
axmu.set_ylabel("$\mu$ (cm$^2/(Vs)$)", fontsize = fontsize)
axmu.legend(fontsize = legend_fontsize)
axmu.tick_params(labelsize = fontsize)
ax2.plot(xTinv, yEa, label = '$E_a(obs)$', color = 'black', linestyle = '-', linewidth = 0.5, marker = 'o')
ax2.plot(xTinv, yEaT23, label = '$E_a$(obs)(log($N$$T^{-2/3}$))', color = 'blue', linestyle = '-', linewidth = 0.5, marker = '^')
ax2.plot(xTinv, yEaT43, label = '$E_a$(obs)(log($N$$T^{-4/3}$))', color = 'red', linestyle = '-', linewidth = 0.5, marker = '+')
ax2.set_xlabel("1000/T (K$^{-1}$)", fontsize = fontsize)
ax2.set_ylabel("$E_a$ (eV)", fontsize = fontsize)
ax2.legend(fontsize = legend_fontsize)
ax2.tick_params(labelsize = fontsize)
ax4.plot(xT, yEa, label = '$E_a(obs)$', color = 'black', linestyle = '-', linewidth = 0.5, marker = 'o')
ax4.plot(xT, yEaT23, label = '$E_a$(obs)(log($N$$T^{-2/3}$))', color = 'blue', linestyle = '-', linewidth = 0.5, marker = '^')
ax4.plot(xT, yEaT43, label = '$E_a$(obs)(log($N$$T^{-4/3}$))', color = 'red', linestyle = '-', linewidth = 0.5, marker = '+')
ax4.set_xlabel("T (K)", fontsize = fontsize)
ax4.set_ylabel("$E_a$ (eV)", fontsize = fontsize)
ax4.legend(fontsize = legend_fontsize)
ax4.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
app.terminate("", usage = usage)
[ドキュメント]
def main(app, cparams, dos, mobility):
"""
概要: アプリケーションのメイン実行関数です。
詳細説明:
初期化されたアプリケーション、パラメータ、DOS、移動度オブジェクトを受け取り、
コマンドライン引数と設定ファイルに基づいて設定を更新します。
初期パラメータを表示し、バンド構造や不純物準位に関する詳細な情報も出力します。
その後、`cparams.mode` で指定されたモード（シミュレーションまたはフィッティング）
に応じて、適切な解析関数 (`sim`, `fit_RH`, `fit`) を呼び出します。
:param app: アプリケーションオブジェクト。
:type app: tkApplication
:param cparams: アプリケーションパラメータオブジェクト。
:type cparams: tkParams
:param dos: 状態密度オブジェクト。
:type dos: tkDOS
:param mobility: 移動度オブジェクト。
:type mobility: tkMobility
:returns: None
:rtype: None
"""
updatevars(app, cparams, dos, mu)
print("")
print("===============================================================================")
print(" Calculate T dependences of semiconductor statistics using VASP results")
print("===============================================================================")
print("")
print("infile : {}".format(cparams.infile))
print("parameter file : {}".format(cparams.parameterfile))
print("parameter backup file: {}".format(cparams.parameterbkfile))
print("mode: ", cparams.mode)
print("")
print("")
print("Read [{}]".format(cparams.parameterfile))
cparams.read_parameters(cparams.parameterfile)
dos.read_parameters(cparams.parameterfile, dos.optid, mobility)
dos.set_misc_parameters(cparams.get_dict())
# check args again to use given parameters
updatevars(app, cparams, dos, mobility)
print("")
print("Initial parameters")
print_parameters(dos.varname, dos.varunit, dos.parameter_list(), dos.optid)
dos.Eg = dos.EC - dos.EV
dos.NC = dos.meff2NC_FEA(dos.meeff, cparams.T0)
dos.NV = dos.meff2NC_FEA(dos.mheff, cparams.T0)
dos.DC0 = dos.meff2DC0_FEA(dos.meeff, cparams.T0)
dos.DV0 = dos.meff2DC0_FEA(dos.mheff, cparams.T0)
print("")
print("Eg: {} eV".format(dos.Eg))
print("Valence band:")
print(" EV : {} eV".format(dos.EV))
print(" mh*: {} me".format(dos.mheff))
print(" NV : {} me at {} K".format(dos.NV, cparams.T0))
print(" DV0: {} me at {} K".format(dos.DV0, cparams.T0))
print("Conduction band:")
print(" EC : {} eV".format(dos.EC))
print(" me*: {} me".format(dos.meeff))
print(" NC : {} me at {} K".format(dos.NC, cparams.T0))
print(" DC0: {} me at {} K".format(dos.DC0, cparams.T0))
print("")
print("Acceptor states:")
print(" {} cm^-3 at {} eV".format(dos.NA, dos.EA))
print("Donoar states:")
print(" {} cm^-3 at {} eV".format(dos.ND, dos.ED))
print("")
print("Scattering factor: r = {}".format(mobility.rfac))
dos.EFmin = dos.EV + cparams.dEFmin
dos.EFmax = dos.EC + cparams.dEFmax
print("d=", dos.EFmin, dos.EFmax,dos.nEF)
dos.EFstep = (dos.EFmax - dos.EFmin) / (dos.nEF - 1)
dos.Emin = dos.EV + cparams.dEFmin
dos.Emax = dos.EC + cparams.dEFmax
dos.nE = int((dos.Emax - dos.Emin) / dos.Estep + 1.00001)
dos.Elist = [dos.Emin + i * dos.Estep for i in range(dos.nE)]
dos.DOSlist = [dos.DOS(dos.Elist[i]) for i in range(dos.nE)]
print("")
print("DOS range: {:8.4g} - {:8.4g}, {:8.4g} steps, nE={}".format(dos.Emin, dos.Emax, dos.Estep, dos.nE))
print("")
print("Integration configuration")
print(" Integration E range: {} - {} eV, or EF+-{}*kBT eV".format(dos.Einteg0, dos.Einteg1, dos.nrange))
if cparams.mode == 'sim':
sim(app, cparams, dos, mobility)
elif cparams.mode == 'fit_RH':
fit_RH(app, cparams, dos, mobility)
elif 'fit' in cparams.mode:
fit(app, cparams, dos, mobility)
else:
app.terminate("Error in main: Invalide cparams.mode [{}]".format(cparams.mode), usage = usage)
if __name__ == '__main__':
app, cparams, dos, mu = initialize()
main(app, cparams, dos, mu)