Seebeck2.py ダウンロード/コピー

Seebeck2.py をダウンロード

Seebeck2.py
Seebeck2.py
   1"""
   2概要:
   3    Seebeck係数と温度から有効質量を推定するスクリプトです。
   4詳細説明:
   5    G.J. Snyder, A. Pereyra and R. Gurunathan, Adv. Funct. Mater. 32, 2112772 (2022)
   6    の論文に記載された手法に基づいて、Seebeck係数から有効質量を推定します。
   7    本スクリプトは、フェルミ積分やバンド構造の計算を含む複数のモードをサポートし、
   8    材料の熱電特性、特にゼーベック係数、電気伝導率、移動度、ローレンツ数、力率、ZT値を計算、
   9    シミュレーション、および実験データへのフィッティングを行います。
  10    Hall係数やHall因子、Hall移動度の計算も行います。
  11主な機能:
  12    - 基本関数 (フェルミ・ディラック分布、フェルミ積分) の計算とプロット。
  13    - 熱電特性 (Seebeck係数、電気伝導率、移動度、熱伝導率、力率、ZT値) の計算とプロット。
  14    - Hall特性 (Hall係数、Hall因子、Hall移動度) の計算とプロット。
  15    - 実験データへのパラメータフィッティング (有効質量、平均自由行程因子)。
  16    - 温度依存性のシミュレーション。
  17    - ローレンツ数と電子熱伝導率の計算。
  18    - Hall因子とドリフト移動度、キャリア濃度の計算。
  19関連リンク:
  20    Seebeck2_usage
  21"""
  22import os
  23import sys
  24import openpyxl
  25from math import exp, sqrt, log, gamma
  26import numpy as np
  27#from numpy import arange
  28from scipy import optimize          # newton関数はscipy.optimizeモジュールに入っている
  29from scipy.interpolate import interp1d
  30from matplotlib import pyplot as plt
  31from scipy import optimize          # newton関数はscipy.optimizeモジュールに入っている
  32from scipy.optimize import minimize
  33
  34
  35from tklib.tkutils import terminate, getarg, getintarg, getfloatarg, safe_getelement, validate_error
  36from tklib.tkutils import pint, pfloat, sort_lists
  37from tklib.tkinifile import tkIniFile
  38from tklib.tkvariousdata import tkVariousData
  39from tklib.tksci.tksci import pi, h, hbar,c, e, kB, me
  40from tklib.tksci.tkmatrix import make_matrix1, make_matrix2
  41from tklib.tkapplication import tkApplication
  42
  43from tklib.tktransport.tkTransport import read_datafile, FermiIntegral_x2
  44from tklib.tktransport.tkTransport import FermiIntegral_fast as FermiIntegral
  45
  46try:
  47    import tklib.tktransport.tkfermi_integral
  48    from tklib.tktransport.tkfermi_integral import FermiIntegral_fast as FermiIntegral_fast_new, fermi_integral_type
  49except:
  50    FermiIntegral_fast_new = None
  51
  52from tklib.tktransport.tkTransport import integrate_Simpson_list, fe, fh, meff2NC_FEA, meff2DC0_FEA, BMShift_FEA
  53from tklib.tktransport.tkWeightedMobility import weighted_mobility, weighted_mobility_new, weighted_mobility_exact
  54from tklib.tktransport.tkTransport import LorentzNumber_FEA
  55from tklib.tktransport.tkDOS_FEA import integrate_Simpson, integrate_Simpson_list, tkDOS
  56from tklib.tktransport.tkmobility_tau import tkMobility, split_optstr
  57
  58
  59#================================
  60# Global variables
  61#================================
  62Debug = 0
  63
  64if FermiIntegral_fast_new:
  65    fermi_integral_type = tkfermi_integral.fermi_integral_type
  66else:
  67    fermi_integral_type = "build-in"
  68
  69L_FEA = LorentzNumber_FEA
  70
  71#mode: 'basic', 'prop', 'Hall', 'init', 'sim', 'fit', 'T', 'calL', 'calF'
  72mode = 'sim'
  73
  74infile      = 'SnSeTe-S-Hall.xlsx'
  75T_label     = r'^T[\s|\(|\[|$]'
  76sigma_label = r'^sigma[$|\S]?.*$'
  77n_label     = r'^N[$|\S]?.*$'
  78mu_label    = r'^mu[$|\s|\(\[]?.*$'
  79S_label     = r'^S[$|\s|\(\[]?.*$'
  80
  81outxlsFjfile   = "Fj.xlsx"
  82outxlsfFDfile  = "fFD.xlsx"
  83outxlsPropfile = "properties.xlsx"
  84outxlsSLfile   = "S-L.xlsx"
  85outcsvFitfile  = ""
  86outcsvLfile    = ""
  87
  88parameterfile   = None
  89parameterbkfile = None
  90outxlsfile      = None
  91
  92# Calculation ranges for mode = 'basic'
  93# x = E / kB / T
  94xmin_basic  = -20.0
  95xmax_basic  =  20.0
  96nx_basic    = 401
  97xstep_basic = (xmax_basic - xmin_basic) / (nx_basic - 1)
  98
  99# Calculation ranges for mode = 'prop'
 100# x = E / kB / T
 101xmin  = -30.0
 102xmax  = 250.0
 103nx    = 281
 104xstep = (xmax - xmin) / (nx - 1)
 105
 106# Seebeck coefficient
 107Smin  = 0.0
 108Smax  = 1.0e-3  # V/K
 109nS    = 101
 110Sstep = (Smax - Smin) / (nS - 1)
 111
 112# N range
 113Nmin = 1.0e15  # cm-3
 114Nmax = 1.0e22  # cm-3
 115nN   = 101
 116
 117# T range
 118T0 = 300.0   # K
 119Tmin = 300.0
 120Tmax = 800.0
 121nT   = 6
 122
 123# read from DOSCAR
 124dos = tkDOS()
 125dos.meeff = 0.3
 126#dos.mheff = 1.0
 127#dos.EV = 0.0
 128dos.EC = 0.0
 129#dos.EA = 0.05   # Acceptor level, eV
 130#dos.NA = 0.0e17
 131#dos.ED = 1.05   # Donor level
 132#dos.ND = 0.0e17
 133#dos.EF0 = 0.0   # initial EF to find EF
 134
 135
 136#移動度パラメータ
 137mobility = tkMobility()
 138# For debug purpose. 1 will use 3-terms polynomial for the inverse of mobility
 139mobility.debug      = 0
 140mobility.use_simple = 0
 141#mobility.use_simple = 1
 142# charge, effective mass, scattering factor
 143mobility.charge = 1.0            # in e, 1.0 for hole, -1.0 for electron
 144mobility.meff   = dos.meeff      # in me
 145mobility.rfac   = 0.5            # tau = (meff/2)^0.5 * l0(T) * E^(r-0.5)
 146mobility.l0     = 1.0e-8         # m
 147
 148# Lattice thermal conductivity
 149klatt  = 5.0            # W/m/K
 150
 151# EF range for mode = 'EF'
 152#dEFmin    = -1.0  # measured from EV, eV
 153#dEFmax    =  1.0  # measured from EC, eV
 154#dos.nEF   = 50
 155#dos.Estep = 0.01    # Integration step, eV
 156
 157epsEF = 1.0e-5
 158
 159label_sample = None
 160xsample      = None
 161label_S      = None
 162yS           = None
 163label_sigma = None
 164ysigma       = None
 165label_N      = None
 166yN           = None
 167label_mu     = None
 168ymu          = None
 169
 170ysigmaini = None
 171ymuini    = None
 172ySini     = None
 173
 174# 最適化パラメータの初期値
 175varname = ["meff", "l0"]
 176varunit = [  "me",  "m"]
 177ai0     = []
 178optid   = [     1,    1]
 179
 180
 181app = None
 182
 183
 184#=============================================
 185# scipy.optimize.minimizeで使うアルゴリズム
 186#=============================================
 187#nelder-mead Downhill simplex
 188#powell Modified Powell
 189#cg conjugate gradient (Polak-Ribiere method)
 190#bfgs BFGS法
 191#newton-cg Newton-CG
 192#trust-ncg 信頼領域 Newton-CG 法
 193#dogleg 信頼領域 dog-leg 法
 194#L-BFGS-B’ (see here)
 195#TNC’ (see here)
 196#COBYLA’ (see here)
 197#SLSQP’ (see here)
 198#trust-constr’(see here)
 199#dogleg’ (see here)
 200#trust-exact’ (see here)
 201#trust-krylov’ (see here)
 202method = "nelder-mead"
 203#method = 'cg'
 204#method = 'powell'
 205#method = 'bfgs'
 206
 207maxiter = 1000
 208tol    = 1.0e-4
 209h_diff = 1.0e-3
 210outputinterval = 1
 211
 212Nmin_fit     = '*'
 213Nmax_fit     = '*'
 214sigmamin_fit = '*'
 215sigmamax_fit = '*'
 216
 217#=============================
 218# Graph configuration
 219#=============================
 220fig = None
 221figsize       = (12, 8)
 222figsize_sim   = (10, 6)
 223figsize_small = (8, 6)
 224figsize_FD = (8, 6)
 225fontsize = 18
 226legend_fontsize = 8
 227graphupdateinterval = 10
 228
 229
 230#=============================
 231# Treat argments
 232#=============================
 233def parameter_list():
 234    """
 235    概要:
 236        現在の移動度パラメータのリストを返します。
 237    戻り値:
 238        :returns: 有効質量と平均自由行程プレファクタのタプル。
 239        :rtype: tuple
 240    """
 241    return mobility.meff, mobility.l0
 242
 243def set_parameters(ai):
 244    """
 245    概要:
 246        引数として与えられた値に基づいて、有効質量と平均自由行程のパラメータを設定します。
 247    詳細説明:
 248        dos.meeff、dos.NC、dos.DC0、mobility.meff、mobility.l0 を更新し、
 249        さらに散乱パラメータを再設定します。
 250    引数:
 251        :param ai: 設定するパラメータのリスト。最初の要素は有効質量、2番目の要素は平均自由行程プレファクタです。
 252        :type ai: list
 253    """
 254    global dos, mobility
 255
 256    dos.meeff     = ai[0]
 257    dos.NC        = meff2NC_FEA(dos.meeff, T0)
 258    dos.DC0       = meff2DC0_FEA(dos.meeff, T0)
 259    mobility.meff = ai[0]
 260    if len(ai) > 1:
 261#        print("ai1=", ai[1])
 262        mobility.l0   = ai[1]
 263        mobility.set_scattering_parameters()
 264#        print("l0=", mobility.l0)
 265
 266def read_parameters(path):
 267    """
 268    概要:
 269        指定されたINIファイルからパラメータを読み込み、グローバル変数に設定します。
 270    詳細説明:
 271        ファイルが存在しない場合は処理をスキップします。
 272        主に有効質量、平均自由行程プレファクタ、電荷、散乱因子、EC、T0などの値を読み込みます。
 273    引数:
 274        :param path: パラメータファイルへのパス。
 275        :type path: str
 276    """
 277    global ai0, optid, T0
 278    
 279    ini = tkIniFile()
 280    inf = ini.ReadAll(path, AddSection = 0)
 281    if inf is None:
 282        return
 283
 284    ai0 = list(ai0)
 285    keylist = ["meff", "l0"]
 286    for i in range(2):
 287        key = keylist[i]
 288        str = inf.get(key, None)
 289        val, id = split_optstr(str)
 290        if val is not None:
 291            ai0[i] = val
 292            optid[i] = id
 293
 294    mobility.meff = pfloat(ai0[0])
 295    dos.meeff     = pfloat(ai0[0])
 296    mobility.l0   = pfloat(ai0[1])
 297    
 298    mobility.charge = pfloat(safe_getelement(inf, "charge", mobility.charge))
 299    mobility.rfac   = pfloat(safe_getelement(inf, "rfac", mobility.rfac))
 300    dos.EC          = pfloat(safe_getelement(inf, "EC", dos.EC))
 301    
 302    T0              = pfloat(safe_getelement(inf, "T0", T0))
 303
 304def save_parameters(path, ai, args):
 305    """
 306    概要:
 307        パラメータと引数を指定されたINIファイルに保存します。
 308    引数:
 309        :param path: 保存先のINIファイルへのパス。
 310        :type path: str
 311        :param ai: 保存する最適化パラメータのリスト。
 312        :type ai: list
 313        :param args: 追加で保存する引数の辞書。
 314        :type args: dict
 315    """
 316    ini = tkIniFile(path)
 317    for key in args.keys():
 318        ini.WriteString('Preferences', key, args[key])
 319
 320    for i in range(len(ai)):
 321        ini.WriteString('Parameters', varname[i], "{}:{}".format(ai[i], optid[i]))
 322
 323def save_parameterfile(ai = ai0, S2 = ''):
 324    """
 325    概要:
 326        現在のパラメータとS2値をパラメータファイルに保存します。
 327    引数:
 328        :param ai: 保存する最適化パラメータのリスト。デフォルトはai0です。
 329        :type ai: list
 330        :param S2: 最小二乗法によるフィッティング結果のS2値。デフォルトは空文字列です。
 331        :type S2: str
 332    """
 333    save_parameters(parameterfile, ai, 
 334            {"infile": infile, "T0": T0, 
 335             "charge": mobility.charge, "rfac": mobility.rfac, 
 336             "EC": dos.EC, "NC": dos.NC, "DC0": dos.DC0, 
 337             "S2": S2})
 338
 339def print_parameters(ai = None):
 340    """
 341    概要:
 342        現在のパラメータの値をコンソールに出力します。
 343    引数:
 344        :param ai: 出力するパラメータのリスト。Noneの場合、グローバル変数ai0が使用されます。
 345        :type ai: list or None
 346    """
 347    if ai is None:
 348        ai = ai0
 349    for i in range(len(varname)):
 350        print("  {:10s}: {:14.8g} {:6} optid={}".format(varname[i], ai[i], varunit[i], optid[i]))
 351
 352def usage(app):
 353    """
 354    概要:
 355        スクリプトのコマンドライン引数の使用方法を表示します。
 356    引数:
 357        :param app: アプリケーションオブジェクト。
 358        :type app: tklib.tkapplication.tkApplication
 359    """
 360    print("")
 361    print("Usage: Variables in () are optional")
 362    print(" (i) python {} basic".format(sys.argv[0]))
 363    print("     Plot basic functions")
 364    print(" (ii) python {} prop T meff r l0 kappa_latt".format(sys.argv[0]))
 365    print("          meff in me0, l0 in m, kappa_latt in W/m/K")
 366    print("     Plot Seebeck related properties")
 367    print("     ex: python {} basic {} {} {} {} {}".format(sys.argv[0], T0, dos.meeff, mobility.rfac, mobility.l0, klatt))
 368    print(" (ii') python {} Hall T meff r l0".format(sys.argv[0]))
 369    print("     Plot Hall related properties")
 370    print("     ex: python {} Hall {} {} {} {}".format(sys.argv[0], T0, dos.meeff, mobility.rfac, mobility.l0))
 371    print(" (iii) python {} init infile".format(sys.argv[0]))
 372    print("     Create .in file")
 373    print(" (iv) python {} sim infile T meff r l0".format(sys.argv[0]))
 374    print("     Calculate weighted mobility etc from input file")
 375    print("     Plot Jonker / Pisarenko plots")
 376    print("     ex: python {} sim {} {} {} {} {}".format(sys.argv[0], infile, T0, dos.meeff, mobility.rfac, mobility.l0))
 377    print(" (v) python {} fit T meff r l0".format(sys.argv[0]))
 378    print("     Fit to Jonker / Pisarenko plots")
 379    print("     ex: python {} fit {} {} {} {} {}".format(sys.argv[0], infile, T0, mobility.meff, mobility.rfac, mobility.l0))
 380    print(" (vi) python {} T infile Tmin Tmax nT".format(sys.argv[0]))
 381    print("     Simulate T dependences")
 382    print("     ex: python {} fit {} {} {} {}".format(sys.argv[0], infile, Tmin, Tmax, nT))
 383    print(" (vii) python {} calL infile meff r l0 ".format(sys.argv[0]))
 384    print("     Calculate L and kappa,e from T, Ne and sigma")
 385    print("     ex: python {} calL {} {} {} {}".format(sys.argv[0], infile, dos.meeff, mobility.rfac, mobility.l0))
 386    print(" (viii) python {} calF infile meff r l0 T_label n_label mu_label sigma_label".format(sys.argv[0]))
 387    print("     Calculate FHall from T, Ne and mu")
 388    print("     ex: python {} calF {} {} {} {} {} {} {}".format(sys.argv[0], infile, dos.meeff, mobility.rfac, mobility.l0, 
 389                            T_label, n_label, mu_label, sigma_label))
 390
 391def updatevars():
 392    """
 393    概要:
 394        コマンドライン引数に基づいて、グローバル変数を更新します。
 395    詳細説明:
 396        スクリプトの実行モード (basic, prop, Hall, init, sim, fit, T, calL, calF, help, usage) に応じて、
 397        入力ファイル、ラベル、温度、有効質量、散乱因子、平均自由行程プレファクタなどのパラメータを設定します。
 398        無効なモードが指定された場合はエラーで終了します。
 399    """
 400    global mode, infile, outxlsxfile, parameterfile, parameterbkfile
 401    global T_label, n_label, mu_label, sigma_label, S_label
 402    global mobility, dos, klatt
 403    global T0, Tmin, Tmax, nT
 404    global ai0, optid
 405    global method, tol, maxiter, h_diff
 406    global Nmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fit
 407
 408    argv = sys.argv
 409#    if len(argv) == 1:
 410#        usage()
 411#        exit()
 412
 413    mode   = getarg( 1, mode)
 414    if mode == 'basic':
 415        pass
 416    elif mode == 'prop':
 417        T0              = getfloatarg(2, T0)
 418        dos.meeff       = getfloatarg(3, dos.meeff)
 419        mobility.rfac   = getfloatarg(4, mobility.rfac)
 420        mobility.l0     = getfloatarg(5, mobility.l0)
 421        klatt           = getfloatarg(6, klatt)
 422    elif mode == 'Hall':
 423        T0              = getfloatarg(2, T0)
 424        dos.meeff       = getfloatarg(3, dos.meeff)
 425        mobility.rfac   = getfloatarg(4, mobility.rfac)
 426        mobility.l0     = getfloatarg(5, mobility.l0)
 427    elif mode == 'init':
 428        infile          = getarg     (2, infile)
 429    elif mode == 'sim':
 430        infile          = getarg     ( 2, infile)
 431        T_label         = getarg     ( 3, T_label)
 432        n_label         = getarg     ( 4, n_label)
 433        mu_label        = getarg     ( 5, mu_label)
 434        sigma_label     = getarg     ( 6, sigma_label)
 435        S_label         = getarg     ( 7, S_label)
 436        T0              = getfloatarg( 8, T0)
 437        dos.meeff       = getfloatarg( 9, dos.meeff)
 438        mobility.rfac   = getfloatarg(10, mobility.rfac)
 439        mobility.l0     = getfloatarg(11, mobility.l0)
 440    elif mode == 'fit':
 441        infile          = getarg     (2, infile)
 442        T_label         = getarg     (3, T_label)
 443        n_label         = getarg     (4, n_label)
 444        mu_label        = getarg     (5, mu_label)
 445        sigma_label     = getarg     (6, sigma_label)
 446        S_label         = getarg     (7, S_label)
 447        T0              = getfloatarg(8, T0)
 448        dos.meeff       = getfloatarg(9, dos.meeff)
 449        mobility.rfac   = getfloatarg(10, mobility.rfac)
 450        mobility.l0     = getfloatarg(11, mobility.l0)
 451        method          = getarg     (12, method)
 452        tol             = getfloatarg(13, tol)
 453        maxiter         = getintarg  (14, maxiter)
 454        Nmin_fit        = getarg     (15, Nmin_fit)
 455        Nmax_fit        = getarg     (16, Nmax_fit)
 456        sigmamin_fit    = getarg     (17, sigmamin_fit)
 457        sigmamax_fit    = getarg     (18, sigmamax_fit)
 458    elif mode == 'T':
 459        infile          = getarg     (2, infile)
 460        Tmin            = getfloatarg(3, Tmin)
 461        Tmax            = getfloatarg(4, Tmax)
 462        nT              = getintarg  (5, nT)
 463    elif mode == 'calL' or mode == 'calF':
 464        infile          = getarg     (2, infile)
 465        T_label         = getarg     (3, T_label)
 466        n_label         = getarg     (4, n_label)
 467        mu_label        = getarg     (5, mu_label)
 468        sigma_label     = getarg     (6, sigma_label)
 469        dos.meeff       = getfloatarg(7, dos.meeff)
 470        mobility.rfac   = getfloatarg(8, mobility.rfac)
 471        mobility.l0     = getfloatarg(9, mobility.l0)
 472    elif mode == 'help' or mode == 'usage':
 473        app.terminate("", usage = usage, pause = True)
 474    else:
 475        app.terminate("Error in updatevars(): Invalid mode {}".format(mode), usage = usage, pause = True)
 476
 477    mobility.meff = dos.meeff
 478    mobility.set_scattering_parameters(mobility.l0, mobility.rfac)
 479    
 480    ai0 = parameter_list()
 481
 482
 483def basic():
 484    """
 485    概要:
 486        フェルミ・ディラック分布とフェルミ積分の基本計算を実行し、結果をファイルに出力し、プロットします。
 487    詳細説明:
 488        x軸を E / kB / T として、電子のフェルミ・ディラック分布、正孔のフェルミ・ディラック分布、
 489        -df/dx、およびそれらの近似値を計算します。
 490        また、フェルミ積分 F_r(x) (r=0, 0.5, 1.0, 1.5, 2.0) の値を計算し、理論値 (exp(x) * Gamma(r+1)) と比較します。
 491        計算結果はoutxlsfFDfileとoutxlsFjfileに保存され、matplotlibでプロットされます。
 492    """
 493    global mobility, dos
 494    
 495    print("")
 496    print("mode:", mode)
 497
 498    print("")
 499    print("Fermi-Dirac distribution")
 500    xx         = []
 501    yfFD       = []
 502    yfFDh      = []
 503    ymdfdx     = []
 504    yfFDapprox = []
 505    print("{:8}\t{:12}\t{:12}\t{:12}\t{:12}".format("x=E/kBT", "fFDe(exact)", "fFDe(approx)", "fFDh(exact)", "-df/dx"))
 506    for i in range(nx_basic):
 507        x  = xmin_basic + i * xstep_basic
 508        xx.append(x)
 509        fFD   = 1.0 / (exp(x) + 1.0)
 510        fFDh  = 1.0 / (exp(-x) + 1.0)
 511        mdfdx = fFD * fFDh
 512        fFDa  = 0.5 - 0.25 * x + 1.0 / 48.0 * x**3 - 1.0 / 64.0 * x**5
 513        
 514        yfFD.append(fFD)
 515        yfFDh.append(fFDh)
 516        ymdfdx.append(mdfdx)
 517        yfFDapprox.append(fFDa)
 518
 519        print("{:8.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}".format(x, fFD, fFDa, fFDh, mdfdx))
 520
 521    print("")
 522    print("Fermi integrals")
 523    yF00 = []
 524    yF05 = []
 525    yF10 = []
 526    yF15 = []
 527    yF20 = []
 528    yG00 = []
 529    yG05 = []
 530    yG10 = []
 531    yG15 = []
 532    yG20 = []
 533    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)"))
 534    for i in range(nx_basic):
 535        x  = xx[i]
 536        for j in range(3):
 537            eta = j / 2.0
 538            y0 = FermiIntegral(x, j)
 539            y1 = FermiIntegral_x2(x, j)
 540            validate_error(y0, y1, 2.0e-6, "Error in basic() for F{}: ".format(eta))
 541
 542        yF00.append(FermiIntegral(x, 0.0))
 543        yF05.append(FermiIntegral(x, 0.5))
 544        yF10.append(FermiIntegral(x, 1.0))
 545        yF15.append(FermiIntegral(x, 1.5))
 546        yF20.append(FermiIntegral(x, 2.0))
 547        yG00.append(exp(x) * gamma(1.0))
 548        yG05.append(exp(x) * gamma(1.5))
 549        yG10.append(exp(x) * gamma(2.0))
 550        yG15.append(exp(x) * gamma(2.5))
 551        yG20.append(exp(x) * gamma(3.0))
 552        if i % 10 == 0:
 553            print("{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}\t{:12.4g}"
 554                    .format(x, yF00[i], yF05[i], yF10[i], yF15[i], yF20[i], yG00[i]))
 555
 556    print("")
 557    print("Save data to [{}]".format(outxlsfFDfile))
 558    tkVariousData().to_excel(outxlsfFDfile, ["x=E/kBT", "fFD(exact)", "fFD(approx)", "fFDh(exact)", "-df/dx"], 
 559                            [xx, yfFD, yfFDapprox, yfFDh, ymdfdx])
 560
 561    print("Save data to [{}]".format(outxlsFjfile))
 562    tkVariousData().to_excel(outxlsFjfile, ["x=EF/kBT", "F0", "F1/2", "F1", "F3/2", "F2", "exp(x)*G(1)", 
 563                            "exp(x)*G(1.5)", "exp(x)*G(2)", "exp(x)*G(2.5)", "exp(x)*G(3)"], 
 564                            [xx, yF00, yF05, yF10, yF15, yF20, yG00, yG05, yG10, yG15, yG20])
 565
 566
 567#=============================
 568# グラフの表示
 569#=============================
 570    print("")
 571
 572    fig = plt.figure(figsize = figsize_FD)
 573
 574    ax2   = fig.add_subplot(2, 1, 1)
 575    ax3   = fig.add_subplot(2, 1, 2)
 576
 577    ax2.plot(xx, yfFD,       label = '$f_{FD,e}(exact)$', linestyle = '-',      color = 'black', linewidth = 0.5)
 578    ax2.plot(xx, yfFDh,      label = '$f_{FD,h}(exact)$', linestyle = '-',      color = 'blue', linewidth = 0.5)
 579    ax2.plot(xx, yfFDapprox, label = '$f_{FD}(approx)$',  linestyle = '-',      color = 'green', linewidth = 0.5)
 580    ax2.plot(xx, ymdfdx,     label = '$-df/dx$',          linestyle = 'dashed', color = 'red', linewidth = 0.5)
 581    ax2.set_xlabel("$x=(E-E_F)/k_BT$ (eV)", fontsize = fontsize)
 582    ax2.set_ylabel("$f_{FD}$, $-df/dx$", fontsize = fontsize)
 583    ax2.set_xlim([-10.0, 10.0])
 584    ax2.set_ylim([-0.1, 1.1])
 585    ax2.legend(fontsize = legend_fontsize)
 586    ax2.tick_params(labelsize = fontsize)
 587
 588    ax3.plot(xx, yF00, label = '$F_0$',     linestyle = '-', color = 'black', linewidth = 0.5)
 589    ax3.plot(xx, yF05, label = '$F_{1/2}$', linestyle = '-', color = 'red',   linewidth = 0.5)
 590    ax3.plot(xx, yF10, label = '$F_1$',     linestyle = '-', color = 'blue',  linewidth = 0.5)
 591    ax3.plot(xx, yF15, label = '$F_{3/2}$', linestyle = '-', color = 'purple', linewidth = 0.5)
 592    ax3.plot(xx, yF20, label = '$F_2$',     linestyle = '-', color = 'green', linewidth = 0.5)
 593    ax3.plot(xx, yG00, label = r'$e^x$$\Gamma(1)$',   linestyle = '', marker = 'o', markerfacecolor = 'black',  markersize = 1)
 594    ax3.plot(xx, yG05, label = r'$e^x$$\Gamma(3/2)$', linestyle = '', marker = 'o', markerfacecolor = 'red',    markersize = 1)
 595    ax3.plot(xx, yG10, label = r'$e^x$$\Gamma(2)$',   linestyle = '', marker = 'o', markerfacecolor = 'blue',   markersize = 1)
 596    ax3.plot(xx, yG15, label = r'$e^x$$\Gamma(5/2)$', linestyle = '', marker = 'o', markerfacecolor = 'purple', markersize = 1)
 597    ax3.plot(xx, yG20, label = r'$e^x$$\Gamma(3)$',   linestyle = '', marker = 'o', markerfacecolor = 'green',  markersize = 1)
 598    ax3.set_xlabel("$x=E_F/k_BT$ (eV)", fontsize = fontsize)
 599    ax3.set_ylabel("$F_r$", fontsize = fontsize)
 600    ax3.set_xlim([-10.0, 10.0])
 601    ax3.set_ylim([1.0e-5, 0.1e4])
 602    ax3.set_yscale('log')
 603    ax3.legend(fontsize = legend_fontsize)
 604    ax3.tick_params(labelsize = fontsize)
 605
 606    plt.tight_layout()
 607    plt.pause(0.1)
 608
 609    app.terminate("", pause = True)
 610
 611def Hall():
 612    """
 613    概要:
 614        ホール係数、ホール因子、ホール移動度などのHall特性を計算し、結果をファイルに出力し、プロットします。
 615    詳細説明:
 616        指定された温度T0と有効質量meff、散乱因子rfac、平均自由行程プレファクタl0に基づいて、
 617        フェルミ準位の掃引に対するキャリア濃度、電気伝導率、移動度、 Hall係数、Hall因子、Hall移動度を計算します。
 618        計算結果はoutxlsPropfileに保存され、matplotlibでプロットされます。
 619    """
 620    print("")
 621    print("Carrier:")
 622    print("  q={} e".format(mobility.charge))
 623    print("  meff={}me".format(dos.meeff))
 624    dos.NC  = meff2NC_FEA(dos.meeff, T0)
 625    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
 626    if mobility.charge < 0.0:
 627        print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
 628        print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
 629    else:
 630        print("  NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
 631        print("  DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
 632    print("  scattering factor (non-degenerated)     r={}".format(mobility.rfac))
 633    print("  mean free path prefactor: l0={} m".format(mobility.l0))
 634
 635    print("")
 636    print(f"Calculat Hall coefficient, Hall factor, and Hall mobility at {T0} K")
 637    xx         = []
 638    yEF        = []
 639    yNe        = []
 640    ysigma     = []
 641    ymu        = []
 642    ytau       = []
 643    yRH0       = []
 644    yRH        = []
 645    yFH        = []
 646    yNe_Hall    = []
 647    ymu_Hall   = []
 648    print("{:>8} {:>8} {:>12} {:>12} {:>12} {:>12} {:>12}"
 649            .format("x=EF/kBT", "EF(eV)", "N(cm^-3)", "sigma(S/cm)", "mu(cm2/Vs)", "FH", "RH(m^3/C"))
 650    for i in range(nx):
 651        x  = xmin + i * xstep
 652        EF = x * kB * T0 / e
 653        
 654        inf = dos.cal_Hall_properteis(T0, EF, mobility, validate_error_str = 'Error in properteis()')
 655
 656        sigma   = inf["sigma"]
 657        Ne      = inf["n"]
 658        mu      = inf["mu"]
 659        tau     = inf["<tau>"]
 660        FH      = inf["FH"]
 661        RH      = inf["RH"]
 662        RH0     = inf["RH0"]
 663        n_Hall  = inf["nHall"]
 664        mu_Hall = inf["muHall"]
 665
 666        xx.append(x)
 667        yEF.append(EF)
 668        yNe.append(Ne)
 669        ysigma.append(sigma)
 670        ymu.append(mu)
 671        ytau.append(tau)
 672        yRH0.append(RH0)
 673        yRH.append(RH)
 674        yFH.append(FH)
 675        yNe_Hall.append(n_Hall)
 676        ymu_Hall.append(mu_Hall)
 677
 678        if i % 10 == 0:
 679            print("{:8.3g} {:8.3g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g}"
 680                    .format(x, EF, yNe[i], ysigma[i], ymu[i], yFH[i], yRH[i]))
 681
 682    print("")
 683    print("Save data to [{}]".format(outxlsPropfile))
 684    tkVariousData().to_excel(outxlsPropfile, ["x=EF/kBT", "EF(eV)", "N(cm^-3)", "sigma(S/cm)", "mu(cm2/Vs)", "FH", "RH(m^3/C"], 
 685                             [xx, yEF, yNe, ysigma, ymu, yFH, yRH])
 686
 687
 688#=============================
 689# グラフの表示
 690#=============================
 691    print("")
 692    fig = plt.figure(figsize = figsize)
 693
 694    axtau   = fig.add_subplot(3, 3, 1)
 695    axNe    = fig.add_subplot(3, 3, 2)
 696    axsigma = fig.add_subplot(3, 3, 3)
 697    axmu    = fig.add_subplot(3, 3, 4)
 698    aRH     = fig.add_subplot(3, 3, 5)
 699    aFH     = fig.add_subplot(3, 3, 6)
 700    aNeRH   = fig.add_subplot(3, 3, 7)
 701    aNeFH   = fig.add_subplot(3, 3, 8)
 702
 703    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))
 704
 705    axtau.plot(xx, ytau,    label = '$\\tau(EF)$',  linestyle = '-', color = 'red',  linewidth = 0.5)
 706#    axtau.plot(xx,     ytauavg, label = '<$\\tau$>',    linestyle = '-', color = 'blue', linewidth = 0.5)
 707    axtau.set_xlabel("$x=(E_F-E_C)/k_BT$ (eV)", fontsize = fontsize)
 708    axtau.set_ylabel("$\\tau$ (fs)", fontsize = fontsize)
 709    axtau.set_ylim([0.0, max(ytau) * 1.1])
 710    axtau.legend(fontsize = legend_fontsize, loc = 'best')
 711    axtau.tick_params(labelsize = fontsize)
 712
 713    axNe.plot(xx, yNe,     label = '$N_e$',           linestyle = '-', color = 'red',  linewidth = 1.0)
 714    axNe.plot(xx, yNe_Hall, label = '$N_e$$_{,Hall}$', linestyle = '-', color = 'blue',  linewidth = 1.0)
 715    axNe.set_xlabel("$x=(E_F-E_C)/k_BT$ (eV)", fontsize = fontsize)
 716    axNe.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
 717    axNe.set_yscale('log')
 718    axNe.legend(fontsize = legend_fontsize, loc = 'best')
 719    axNe.tick_params(labelsize = fontsize)
 720
 721#    axsigma.plot(xx, ysigma, label = r'$\sigma$',    linestyle = '-', color = 'red',  linewidth = 1.0)
 722    axsigma.plot(yNe_Hall, ysigma, label = r'$\sigma$',    linestyle = '-', color = 'red',  linewidth = 1.0)
 723#    axsigma.set_xlabel(r"$x=(E_F-E_C)/k_BT$ (eV)", fontsize = fontsize)
 724    axsigma.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
 725    axsigma.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
 726    axsigma.set_xscale('log')
 727    axsigma.set_yscale('log')
 728    axsigma.set_xlim([1.0e10, 1.0e23])
 729    axsigma.legend(fontsize = legend_fontsize, loc = 'best')
 730    axsigma.tick_params(labelsize = fontsize)
 731
 732#    axmu.plot(xx, ymu, label = r'$\mu$',                linestyle = '-', color = 'red',  linewidth = 1.0)
 733    axmu.plot(yNe_Hall, ymu, label = r'$\mu$',                linestyle = '-', color = 'red',  linewidth = 1.0)
 734#    axmu.plot(xx, ymu_Hall, label = r'$\mu$$_{,Hall}$', linestyle = '-', color = 'blue',  linewidth = 1.0)
 735    axmu.plot(yNe_Hall, ymu_Hall, label = r'$\mu$$_{,Hall}$', linestyle = '-', color = 'blue',  linewidth = 1.0)
 736#    axmu.set_xlabel(r"$x=(E_F-E_C)/k_BT$ (eV)", fontsize = fontsize)
 737    axmu.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
 738    axmu.set_ylabel(r"$\mu$ (cm$^2$/V/s)", fontsize = fontsize)
 739    axmu.set_xscale('log')
 740    axmu.set_xlim([1.0e10, 1.0e23])
 741    axmu.legend(fontsize = legend_fontsize, loc = 'best')
 742    axmu.tick_params(labelsize = fontsize)
 743
 744#    aRH.plot(xx, yRH0, label = r'$R_H$$_0$=$1/e/N_e$', linestyle = '-', color = 'red',  linewidth = 1.0)
 745    aRH.plot(yNe_Hall, yRH0, label = r'$R_H$$_0$=$1/e/N_e$', linestyle = '-', color = 'red',  linewidth = 1.0)
 746#    aRH.plot(xx, yRH,  label = r'$R_H$$_{,Hall}$',     linestyle = '-', color = 'blue',  linewidth = 1.0)
 747    aRH.plot(yNe_Hall, yRH,  label = r'$R_H$$_{,Hall}$',     linestyle = '-', color = 'blue',  linewidth = 1.0)
 748#    aRH.set_xlabel(r"$x=(E_F-E_C)/k_BT$ (eV)", fontsize = fontsize)
 749    aRH.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
 750    aRH.set_ylabel(r"$R_H$ (cm$^3$/C)", fontsize = fontsize)
 751    aRH.set_xscale('log')
 752    aRH.set_yscale('log')
 753    aRH.set_xlim([1.0e10, 1.0e23])
 754    aRH.legend(fontsize = legend_fontsize, loc = 'best')
 755    aRH.tick_params(labelsize = fontsize)
 756
 757#    aFH.plot(xx, yFH, label = r'$F_{Hall}$', linestyle = '-', color = 'red',  linewidth = 1.0)
 758    aFH.plot(yNe_Hall, yFH, label = r'$F_{Hall}$', linestyle = '-', color = 'red',  linewidth = 1.0)
 759#    aFH.set_xlabel(r"$x=(E_F-E_C)/k_BT$ (eV)", fontsize = fontsize)
 760    aFH.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
 761    aFH.set_ylabel("$F_H$", fontsize = fontsize)
 762    aFH.set_xscale('log')
 763    aFH.set_xlim([1.0e10, 1.0e23])
 764    aFH.legend(fontsize = legend_fontsize, loc = 'best')
 765    aFH.tick_params(labelsize = fontsize)
 766
 767    aNeRH.plot(yNe_Hall, yRH, label = r'$R_H$$_{,Hall}$', linestyle = '-', color = 'red',  linewidth = 1.0)
 768    aNeRH.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
 769    aNeRH.set_ylabel(r"$R_H$ (cm$^3$/C)", fontsize = fontsize)
 770    aNeRH.set_xscale('log')
 771    aNeRH.set_yscale('log')
 772    aNeRH.set_xlim([1.0e10, 1.0e23])
 773#    aNeRH.set_ylim([min(yRH) * 0.5, min(yRH) * 1.0e3])
 774#    aNeRH.legend(fontsize = legend_fontsize, loc = 'best')
 775    aNeRH.tick_params(labelsize = fontsize)
 776
 777    aNeFH.plot(yNe_Hall, yFH, label = r'$N_e$', linestyle = '-', color = 'red',  linewidth = 1.0)
 778    aNeFH.set_xlabel(r"$N_e$$_{,Hall}$ (cm$^{-3}$)", fontsize = fontsize)
 779    aNeFH.set_ylabel(r"$F_H$ (cm$^{-3}$)", fontsize = fontsize)
 780    aNeFH.set_xscale('log')
 781    aNeFH.set_xlim([1.0e10, 1.0e23])
 782#    aNeFH.legend(fontsize = legend_fontsize, loc = 'best')
 783    aNeFH.tick_params(labelsize = fontsize)
 784
 785    plt.tight_layout()
 786    plt.pause(0.1)
 787
 788    app.terminate("", pause = True)
 789
 790def properties():
 791    """
 792    概要:
 793        Seebeck係数に関連する熱電特性を計算し、結果をファイルに出力し、プロットします。
 794    詳細説明:
 795        指定された温度T0と有効質量meff、散乱因子rfac、平均自由行程プレファクタl0、格子熱伝導率klattに基づいて、
 796        フェルミ準位の掃引に対するキャリア濃度、電気伝導率、移動度、Seebeck係数、電子熱伝導率、
 797        全熱伝導率、ローレンツ数、力率、ZT値を計算します。
 798        また、非縮退および縮退近似でのSeebeck係数も計算し、比較します。
 799        計算結果はoutxlsPropfileに保存され、matplotlibでプロットされます。
 800    """
 801    global mobility, dos
 802    
 803    print("")
 804    print("mode:", mode)
 805    print("T={} K".format(T0))
 806    print("meff:", dos.meeff)
 807    dos.NC  = meff2NC_FEA(dos.meeff, T0)
 808    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
 809    if mobility.charge < 0.0:
 810        print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
 811        print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
 812    else:
 813        print("  NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
 814        print("  DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
 815    print("  scattering factor in tau(E) prop to E^(r-0.5) (non-degenerated) r={}".format(mobility.rfac))
 816    print("Lorentz number")
 817    print("  Free electron model: {:12.8g} Wohm/K^2".format(L_FEA))
 818    print("Mean free path prefactor: l0={} m".format(mobility.l0))
 819    print("Lattice thermal conductivity: klatt={} W/K".format(klatt))
 820
 821    print("")
 822    print("Properties")
 823    xx         = []
 824    xx_tau     = []
 825    yEF        = []
 826    yNe        = []
 827    ysigma     = []
 828    ymu        = []
 829    ytau       = []
 830    ytauavg    = []
 831    yS         = []
 832    ySndeg     = []
 833    ySdeg      = []
 834    ykappa     = []
 835    ykappatot  = []
 836    yL         = []
 837    yLapprox   = []
 838    yPF        = []
 839    yZT        = []
 840    print("{:>8} {:>8} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12}"
 841            .format("x=EF/kBT", "EF(eV)", "N(cm^-3)", "sigma(S/cm)", "mu(cm2/Vs)", "<tau>(fs)",
 842                    "S(uV/K)", "Snon-deg", "Sdeg", "kappa,e(W/m/K)", "kappa,tot", "L(Wohm/K^2)", "Lapprox", "PF(uW/cm/K^2)", "ZT"))
 843    for i in range(nx):
 844        x  = xmin + i * xstep
 845        EF = x * kB * T0 / e
 846        
 847#        print("T=", T0)
 848        sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
 849            = dos.cal_transport_properteis(T0, EF, mobility, klatt, validate_error_str = 'Error in properteis()')
 850
 851# Relaxation time tau(EF)
 852        if EF <= 0.0:
 853            tau = None
 854        else:
 855            tau = mobility.tau(T0, EF) * 1.0e15      # fs, E in eV
 856        if EF > 0.0:
 857            xx_tau.append(x)
 858            ytau.append(tau)
 859
 860        Sndeg = dos.cal_S_nondegenerated_from_Ne(n, mobility.rfac, charge = mobility.charge) * 1.0e6    # uV/K
 861        Sdeg  = dos.cal_S_degenerated_from_Ne(T0, n, mobility.rfac, charge = mobility.charge) * 1.0e6  # uV/K
 862
 863# Approximated Lorentz number
 864        Lapprox = dos.cal_LorentzNumber_from_S_approximate(S)
 865
 866        xx.append(x)
 867        yEF.append(EF)
 868        yNe.append(n)
 869        ysigma.append(sigma)
 870        ymu.append(mu)
 871        ytauavg.append(tau_avg)
 872        yS.append(S)
 873        ySndeg.append(Sndeg)
 874        ySdeg.append(Sdeg)
 875        ykappa.append(kappa)
 876        ykappatot.append(kappa + klatt)
 877        yL.append(L)
 878        yLapprox.append(Lapprox)
 879        yPF.append(PF)
 880        yZT.append(ZT)
 881        if i % 10 == 0:
 882            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}"
 883                    .format(x, EF, n, sigma, mu, tau_avg, S, Sndeg, Sdeg, kappa, kappa_tot, L, Lapprox, PF, ZT))
 884
 885    print("")
 886    print("Save data to [{}]".format(outxlsPropfile))
 887    tkVariousData().to_excel(outxlsPropfile, ["x=EF/kBT", "EF(eV)", "N(cm^-3)", "sigma(S/cm)", "mu(cm2/Vs)", "<tau>(fs)",
 888                    "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"], 
 889                             [xx, yEF, yNe, ysigma, ymu, ytauavg, yS, ySndeg, ySdeg, ykappa, ykappatot, yL, yLapprox, yPF, yZT])
 890
 891
 892#=============================
 893# グラフの表示
 894#=============================
 895    print("")
 896    maxS = max(yS)
 897    minS = min(yS)
 898    if maxS > 0.0:
 899        maxS *= 1.2
 900    else:
 901        maxS = 0.0
 902    if minS > 0.0:
 903        minS = 0.0
 904    else:
 905        minS *= 2.0
 906
 907    fig = plt.figure(figsize = figsize)
 908
 909    axtau   = fig.add_subplot(3, 4, 1)
 910    axNe    = fig.add_subplot(3, 4, 2)
 911    axsigma = fig.add_subplot(3, 4, 3)
 912    axmu    = fig.add_subplot(3, 4, 4)
 913    axS     = fig.add_subplot(3, 4, 5)
 914    axke    = fig.add_subplot(3, 4, 6)
 915    axPF    = fig.add_subplot(3, 4, 7)
 916    axZT    = fig.add_subplot(3, 4, 8)
 917    axL1    = fig.add_subplot(3, 4, 9)
 918    axL2    = fig.add_subplot(3, 4, 10)
 919    axLS    = fig.add_subplot(3, 4, 11)
 920
 921    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))
 922
 923    axtau.plot(xx_tau, ytau,    label = '$\\tau(EF)$',  linestyle = '-', color = 'red',  linewidth = 0.5)
 924    axtau.plot(xx,     ytauavg, label = '<$\\tau$>',    linestyle = '-', color = 'blue', linewidth = 0.5)
 925    axtau.set_xlabel("$E$ (eV)", fontsize = fontsize)
 926    axtau.set_ylabel("$\\tau$ (fs)", fontsize = fontsize)
 927    axtau.set_ylim([0.0, max(ytau) * 1.1])
 928    axtau.legend(fontsize = legend_fontsize, loc = 'best')
 929    axtau.tick_params(labelsize = fontsize)
 930
 931    axNe.plot(xx, yNe, label = '$N_e$',    linestyle = '-', color = 'red',  linewidth = 1.0)
 932    axNe.set_xlabel("$x=(E_F-E_C)/k_BT$ (eV)", fontsize = fontsize)
 933    axNe.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
 934    axNe.set_yscale('log')
 935    axNe.legend(fontsize = legend_fontsize, loc = 'best')
 936    axNe.tick_params(labelsize = fontsize)
 937
 938#    axsigma.plot(xx, ysigma, label = r'$\sigma$',    linestyle = '-', color = 'red',  linewidth = 1.0)
 939    axsigma.plot(yNe, ysigma, label = r'$\sigma$',    linestyle = '-', color = 'red',  linewidth = 1.0)
 940#    axsigma.set_xlabel(r"$x=(E_F-E_C)/k_BT$ (eV)", fontsize = fontsize)
 941    axsigma.set_xlabel(r"$N_e$ (cm$^{-3}$)", fontsize = fontsize)
 942    axsigma.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
 943    axsigma.set_xscale('log')
 944    axsigma.set_yscale('log')
 945    axsigma.legend(fontsize = legend_fontsize, loc = 'best')
 946    axsigma.tick_params(labelsize = fontsize)
 947
 948#    axmu.plot(xx, ymu, label = '$\\mu$',    linestyle = '-', color = 'red',  linewidth = 1.0)
 949    axmu.plot(yNe, ymu, label = '$\\mu$',    linestyle = '-', color = 'red',  linewidth = 1.0)
 950#    axmu.set_xlabel("$x=(E_F-E_C)/k_BT$ (eV)", fontsize = fontsize)
 951    axmu.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
 952    axmu.set_ylabel("$\\mu$ (cm$^2$/V/s)", fontsize = fontsize)
 953    axmu.set_xscale('log')
 954    axmu.legend(fontsize = legend_fontsize, loc = 'best')
 955    axmu.tick_params(labelsize = fontsize)
 956
 957    axke.plot(yNe, ykappa, label = '$\\kappa$$_e$',    linestyle = '-', color = 'red',  linewidth = 1.0)
 958    axke.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
 959#    axke.set_xlabel("$x=(E_F-E_C)/k_BT$ (eV)", fontsize = fontsize)
 960    axke.set_ylabel("$\\kappa$$_e$ (W/m/K)", fontsize = fontsize)
 961    axke.set_xscale('log')
 962    axke.set_yscale('log')
 963    axke.legend(fontsize = legend_fontsize, loc = 'best')
 964    axke.tick_params(labelsize = fontsize)
 965
 966    axS.plot(yNe, yS,     label = '$S$',           linestyle = '-',      color = 'black', linewidth = 1.0)
 967    axS.plot(yNe, ySndeg, label = '$S_{non-deg}$', linestyle = 'dashed', color = 'red',   linewidth = 0.5)
 968    axS.plot(yNe, ySdeg,  label = '$S_{deg}$',     linestyle = 'dashed', color = 'blue',  linewidth = 0.5)
 969    axS.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
 970    axS.set_ylabel(r"$S$ ($\mu$V/K)", fontsize = fontsize)
 971    axS.set_xscale('log')
 972    axS.set_ylim([minS, maxS])
 973    axS.legend(fontsize = legend_fontsize, loc = 'best')
 974    axS.tick_params(labelsize = fontsize)
 975
 976    axPF.plot(yNe, yPF, label = '$PF$',           linestyle = '-',      color = 'black', linewidth = 1.0)
 977    axPF.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
 978    axPF.set_ylabel(r"$PF$ ($\mu$W/cm/K$^2$)", fontsize = fontsize)
 979    axPF.set_xscale('log')
 980    axPF.legend(fontsize = legend_fontsize, loc = 'best')
 981    axPF.tick_params(labelsize = fontsize)
 982
 983    axZT.plot(yNe, yZT, label = '$ZT$',           linestyle = '-',      color = 'black', linewidth = 1.0)
 984    axZT.set_xlabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
 985    axZT.set_ylabel("$ZT$", fontsize = fontsize)
 986    axZT.set_xscale('log')
 987    axZT.legend(fontsize = legend_fontsize, loc = 'best')
 988    axZT.tick_params(labelsize = fontsize)
 989
 990#    axL1.plot(xx, yL, label = 'L(cal)', linestyle = '-', linewidth = 0.5)
 991    axL1.plot(yNe, yL, label = 'L(cal)', linestyle = '-', linewidth = 0.5)
 992    axL1.plot([min(xx), max(xx)], [L_FEA, L_FEA], label = 'L(FEA deg)', linestyle = '-', color = 'red', linewidth = 0.5)
 993#    axL1.set_xlabel("$x=(E_F-E_C)/k_B/T$", fontsize = fontsize)
 994    axL1.set_xlabel(r"$N_e$ (cm$^{-3}$)", fontsize = fontsize)
 995    axL1.set_ylabel(r"$L=\kappa_e/\sigma/T$ (W$\Omega$/K$^2$)", fontsize = fontsize)
 996    axL1.set_ylim([0.0, 3.2e-8])
 997    axL1.legend(fontsize = legend_fontsize, loc = 'best')
 998    axL1.tick_params(labelsize = fontsize)
 999
1000    axL2.plot(yNe, yL, label = 'L(cal)', linestyle = '-', linewidth = 0.5)
1001    axL2.plot([min(yNe), max(yNe)], [L_FEA, L_FEA], label = 'L(FEA deg)', linestyle = '-', color = 'red', linewidth = 0.5)
1002    axL2.set_xlabel(r"$N_e$ (cm$^{-3}$)", fontsize = fontsize)
1003    axL2.set_ylabel(r"$L=\kappa_e/\sigma/T$ (W$\Omega$/K$^2$)", fontsize = fontsize)
1004    axL2.set_xscale('log')
1005    axL2.set_xlim([1.0e10, 1.0e23])
1006    axL2.set_ylim([0.0, 3.2e-8])
1007    axL2.legend(fontsize = legend_fontsize, loc = 'best')
1008    axL2.tick_params(labelsize = fontsize)
1009
1010    axLS.plot(yS, yL,       label = 'L(cal)',         linestyle = '-', color = 'red',  linewidth = 1.0)
1011    axLS.plot(yS, yLapprox, label = 'L(approx, r=0)', linestyle = '-', color = 'blue', linewidth = 0.5)
1012    axLS.plot([min(yS), max(yS)], [L_FEA, L_FEA], label = 'L(FEA deg)', linestyle = '-', color = 'red', linewidth = 0.5)
1013    axLS.set_xlabel(r"$S$ ($\mu$V/K)", fontsize = fontsize)
1014    axLS.set_ylabel(r"$L=\kappa_e/\sigma/T$ (W$\Omega$/K$^2$)", fontsize = fontsize)
1015    axLS.set_ylim([0.0, 3.2e-8])
1016    axLS.legend(fontsize = legend_fontsize, loc = 'best')
1017    axLS.tick_params(labelsize = fontsize)
1018
1019    plt.tight_layout()
1020    plt.pause(0.1)
1021
1022    app.terminate("", pause = True)
1023
1024def recover_parameters(ais, optid, aidef0):
1025    """
1026    概要:
1027        最適化されたパラメータのサブセットから完全なパラメータリストを再構築します。
1028    引数:
1029        :param ais: 最適化されたパラメータのサブセット。
1030        :type ais: list
1031        :param optid: 各パラメータが最適化対象かどうかを示すIDのリスト (1は最適化対象、0は固定)。
1032        :type optid: list
1033        :param aidef0: デフォルトのパラメータ値のリスト。
1034        :type aidef0: list
1035    戻り値:
1036        :returns: 再構築された完全なパラメータリスト。
1037        :rtype: list
1038    """
1039    aidef = list(aidef0).copy()
1040    
1041    ai = []
1042    c = 0
1043    for i in range(len(optid)):
1044        if optid[i] == 1:
1045            ai.append(ais[c])
1046            c += 1
1047        else:
1048            ai.append(aidef[i])
1049
1050    return ai
1051
1052def choose_parameters(ais0, optid):
1053    """
1054    概要:
1055        完全なパラメータリストから最適化対象のパラメータのみを抽出します。
1056    引数:
1057        :param ais0: 完全なパラメータリスト。
1058        :type ais0: list
1059        :param optid: 各パラメータが最適化対象かどうかを示すIDのリスト (1は最適化対象、0は固定)。
1060        :type optid: list
1061    戻り値:
1062        :returns: 最適化対象のパラメータのみを含むリスト。
1063        :rtype: list
1064    """
1065    ais = list(ais0).copy()
1066
1067    ai = []
1068    c = 0
1069    for i in range(len(optid)):
1070        if optid[i] == 1:
1071            ai.append(ais[i])
1072
1073    return ai
1074
1075def calS2_sigma(ai):
1076    """
1077    概要:
1078        与えられたパラメータで計算された電気伝導率と実験データの偏差の二乗和 (S2) を計算します。
1079    詳細説明:
1080        主にmobility.l0のフィッティングに用いられます。有効質量はdos.meeffから設定されます。
1081        入力データのNmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fitに基づいてフィッティング範囲を限定します。
1082        各データ点に対してEFを計算し、そこから電気伝導率sigmafinを計算し、実験値との二乗誤差を評価します。
1083    引数:
1084        :param ai: フィッティングするパラメータのリスト。ここではmobility.l0が対象です。
1085        :type ai: list
1086    戻り値:
1087        :returns: 電気伝導率の二乗和S2。
1088        :rtype: float
1089    """
1090    global mu, dos, ysigma
1091
1092    aiorg = parameter_list()
1093
1094    ainew = [ai[0]]
1095    if ainew[0] <= 0.0:
1096        ainew[0] = 1.0e-4
1097
1098    set_parameters([mobility.meff, ainew[0]])
1099
1100    S2 = 0.0
1101    ndata = len(ysigma)
1102    eps = 1.0e-300
1103    for i in range(ndata):
1104        S      = yS[i]
1105        sigma  = ysigma[i]
1106        n      = yN[i]
1107        mu     = ymu[i]
1108
1109        if Nmin_fit is not None and n < Nmin_fit:
1110            continue
1111        if Nmax_fit is not None and Nmax_fit < n:
1112            continue
1113        if sigmamin_fit is not None and sigma < sigmamin_fit:
1114            continue
1115        if sigmamax_fit is not None and sigmamax_fit < sigma:
1116            continue
1117
1118        EFfin, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
1119        if not flag:
1120            terminate("Error in callback(): EF calculation did not converge. diffEF={}  eps={} in {} iteration"
1121                    .format(diffEF, epsEF, maxiter))
1122
1123        sigmafin, nfin, mufin, tau_avg, Sfin, kappa, kappa_tot, L, PF, ZT, inf \
1124            = dos.cal_transport_properteis(T0, EFfin, mobility, klatt, validate_error_str = 'Error in fit()')
1125#        print("m,l0=", mobility.meff, mobility.l0)
1126#        print("  sigma=", i, ysigma[i], sigmafin)
1127
1128        dlogsigma = ysigma[i] - sigmafin
1129#        dlogsigma = log(ysigma[i]+eps) - log(sigmafin+eps)
1130        S2 += dlogsigma * dlogsigma
1131
1132    S2 /= (ndata - 1)
1133
1134    set_parameters(aiorg)
1135#    print("S2=", S2)
1136    
1137    return S2
1138
1139def calS2_S(ai):
1140    """
1141    概要:
1142        与えられたパラメータで計算されたSeebeck係数と実験データの偏差の二乗和 (S2) を計算します。
1143    詳細説明:
1144        主に有効質量meffのフィッティングに用いられます。
1145        入力データのNmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fitに基づいてフィッティング範囲を限定します。
1146        各データ点に対してEFを計算し、そこからSeebeck係数Sfinを計算し、実験値との二乗誤差を評価します。
1147    引数:
1148        :param ai: フィッティングするパラメータのリスト。ここではdos.meeffが対象です。
1149        :type ai: list
1150    戻り値:
1151        :returns: Seebeck係数の二乗和S2。
1152        :rtype: float
1153    """
1154    global mu, dos, ysigma
1155    global method, tol, maxiter, h_diff
1156    global Nmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fit
1157
1158    aiorg = parameter_list()
1159
1160    ainew = [ai[0]]
1161    if ainew[0] <= 0.0:
1162        ainew[0] = 1.0e-4
1163
1164    set_parameters(ainew)
1165#    print_parameters()
1166#    print("m=", dos.meeff, dos.NC, dos.DC0)
1167
1168    S2 = 0.0
1169    ndata = len(ysigma)
1170    eps = 1.0e-300
1171    for i in range(ndata):
1172        S      = yS[i]
1173        sigma  = ysigma[i]
1174        n      = yN[i]
1175        mu     = ymu[i]
1176
1177        if Nmin_fit is not None and n < Nmin_fit:
1178            continue
1179        if Nmax_fit is not None and Nmax_fit < n:
1180            continue
1181        if sigmamin_fit is not None and sigma < sigmamin_fit:
1182            continue
1183        if sigmamax_fit is not None and sigmamax_fit < sigma:
1184            continue
1185
1186        EFfin, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
1187        if not flag:
1188            app.terminate("Error in callback(): EF calculation did not converge. diffEF={}  eps={} in {} iteration"
1189                    .format(diffEF, epsEF, maxiter), usage = usage, pause = True)
1190
1191        sigmafin, nfin, mufin, tau_avg, Sfin, kappa, kappa_tot, L, PF, ZT, inf \
1192            = dos.cal_transport_properteis(T0, EFfin, mobility, klatt, validate_error_str = 'Error in fit()')
1193        
1194        dlogS = yS[i] - Sfin
1195        """
1196        if yS[i] <= 0.0:
1197            dlogS = log(-yS[i]+eps) - log(-Sfin+eps)
1198        else:
1199            dlogS = log(yS[i]+eps) - log(Sfin+eps)
1200        """
1201        S2 += dlogS * dlogS
1202
1203    S2 /= (ndata - 1)
1204
1205    set_parameters(aiorg)
1206    
1207    return S2
1208
1209# First derivatives to be used e.g. for cg method
1210# Approximate by forward difference method with the delta h = h_diff
1211def diff1_sigma(ai):
1212    """
1213    概要:
1214        電気伝導率の二乗和S2に対するパラメータの1階微分を数値的に計算します (forward difference)。
1215    詳細説明:
1216        最適化アルゴリズム (e.g., 共役勾配法) で使用される勾配の計算に使用されます。
1217        計算されたS2_sigmaの各パラメータに対する微分値を返します。
1218    引数:
1219        :param ai: 評価するパラメータのリスト。
1220        :type ai: list
1221    戻り値:
1222        :returns: S2_sigmaに対する各パラメータの1階微分の配列。
1223        :rtype: numpy.ndarray
1224    """
1225    n = len(ai)
1226    f0 = calS2_sigma(ai)
1227    df = np.empty(n)
1228    for i in range(n):
1229        aii = ai
1230        aii[i] = ai[i] + h_diff
1231        df[i] = (calS2_sigma(aii) - f0) / h_diff
1232    return df
1233
1234def diff1_S(ai):
1235    """
1236    概要:
1237        Seebeck係数の二乗和S2に対するパラメータの1階微分を数値的に計算します (forward difference)。
1238    詳細説明:
1239        最適化アルゴリズム (e.g., 共役勾配法) で使用される勾配の計算に使用されます。
1240        計算されたS2_Sの各パラメータに対する微分値を返します。
1241    引数:
1242        :param ai: 評価するパラメータのリスト。
1243        :type ai: list
1244    戻り値:
1245        :returns: S2_Sに対する各パラメータの1階微分の配列。
1246        :rtype: numpy.ndarray
1247    """
1248    n = len(ai)
1249    f0 = calS2_S(ai)
1250    df = np.empty(n)
1251    for i in range(n):
1252        aii = ai
1253        aii[i] = ai[i] + h_diff
1254        df[i] = (calS2_S(aii) - f0) / h_diff
1255    return df
1256
1257def plot(fig, yn, ysigma, ymu, yS, ysigmaini, ymuini, ySini, 
1258                  ysigmafin = None, ymufin = None, ySfin = None):
1259    """
1260    概要:
1261        PisarenkoプロットとJonkerプロットを表示します。
1262    詳細説明:
1263        実験データ、初期計算値、および最終計算値をそれぞれN vs Sとsigma vs Sの形でプロットします。
1264        フィッティングの進行状況を視覚的に確認するために使用されます。
1265    引数:
1266        :param fig: matplotlibのFigureオブジェクト。
1267        :type fig: matplotlib.figure.Figure
1268        :param yn: キャリア濃度 N のリスト。
1269        :type yn: list
1270        :param ysigma: 電気伝導率 sigma のリスト。
1271        :type ysigma: list
1272        :param ymu: 移動度 mu のリスト。
1273        :type ymu: list
1274        :param yS: Seebeck係数 S のリスト。
1275        :type yS: list
1276        :param ysigmaini: 初期計算された電気伝導率 sigma のリスト。
1277        :type ysigmaini: list
1278        :param ymuini: 初期計算された移動度 mu のリスト。
1279        :type ymuini: list
1280        :param ySini: 初期計算されたSeebeck係数 S のリスト。
1281        :type ySini: list
1282        :param ysigmafin: 最終計算された電気伝導率 sigma のリスト (任意)。
1283        :type ysigmafin: list or None
1284        :param ymufin: 最終計算された移動度 mu のリスト (任意)。
1285        :type ymufin: list or None
1286        :param ySfin: 最終計算されたSeebeck係数 S のリスト (任意)。
1287        :type ySfin: list or None
1288    """
1289    global plt
1290    global graphupdateinterval
1291
1292    plt.clf()
1293
1294    ax1  = fig.add_subplot(1, 2, 1)
1295    ax2  = fig.add_subplot(1, 2, 2)
1296    plt.title(infile, fontsize = fontsize)
1297
1298    Srange = [min(yS) * 0.9, max(yS) * 1.1]
1299    yn2, ySi2 = sort_lists([yn, ySini])
1300    ins1 = ax1.plot(yn, yS,       label = 'Pisarenko(obs)', linestyle = 'none', marker = 'o', markerfacecolor = 'blue', markeredgecolor = 'blue')
1301    ins2 = ax1.plot(yn2, ySi2,    label = 'Pisarenko(ini)', linestyle = '-',    color = 'blue', linewidth = 0.5)
1302#    ins2 = ax1.plot(yn, ySini,    label = 'Pisarenko(ini)', linestyle = '-',    color = 'blue', linewidth = 0.5)
1303    if ysigmafin:
1304        yn2, ySf2 = sort_lists([yn, ySini])
1305        ins3 = ax1.plot(yn2, ySf2,    label = 'Pisarenko(fin)', linestyle = '-',      color = 'red')
1306#        ins3 = ax1.plot(yn, ySfin,    label = 'Pisarenko(fin)', linestyle = '-',      color = 'red')
1307    ax1.set_xscale('log')
1308    ax1.set_ylim([0.0, Srange[1]])
1309    ax1.set_xlabel(r'$N_e$ (cm$^{-3}$)', fontsize = fontsize)
1310    ax1.set_ylabel(r'S ($\mu$V/K)', fontsize = fontsize)
1311    ax1.legend(fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
1312    ax1.tick_params(labelsize = fontsize)
1313
1314    ysi, ySi = sort_lists([ysigmaini, ySini])
1315    ins1 = ax2.plot(ysigma,    yS,       label = 'Jonker(obs)', linestyle = 'none', marker = 'o', markerfacecolor = 'blue', markeredgecolor = 'blue')
1316    ins2 = ax2.plot(ysi, ySi,    label = 'Jonker(ini)', linestyle = '-',    color = 'blue', linewidth = 0.5)
1317#    ins2 = ax2.plot(ysigmaini, ySini,    label = 'Jonker(ini)', linestyle = '-',    color = 'blue', linewidth = 0.5)
1318    if ysigmafin:
1319        ysf, ySf = sort_lists([ysigmafin, ySfin])
1320        ins3 = ax2.plot(ysf, ySf,    label = 'Jonker(fin)', linestyle = '-',      color = 'red')
1321#        ins3 = ax2.plot(ysigmafin, ySfin,    label = 'Jonker(fin)', linestyle = '-',      color = 'red')
1322    ax2.set_xscale('log')
1323    ax2.set_ylim([0.0, Srange[1]])
1324    ax2.set_xlabel(r'$\sigma$ (S/cm)', fontsize = fontsize)
1325    ax2.set_ylabel(r'S ($\mu$V/K)', fontsize = fontsize)
1326    ax2.legend(fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
1327    ax2.tick_params(labelsize = fontsize)
1328
1329    plt.tight_layout()
1330    plt.pause(0.01)
1331
1332
1333# Callback function for scipy.optimize.minimize()
1334# Print variables every iteration, and update graph for every graphupdateinterval iterationsおｔ
1335iter = 0
1336def callback(fig, S2):
1337    """
1338    概要:
1339        scipy.optimize.minimizeのコールバック関数として、最適化の進行状況を出力し、グラフを更新します。
1340    詳細説明:
1341        指定された出力間隔でパラメータとS2値をコンソールに表示し、
1342        指定されたグラフ更新間隔でJonkerプロットとPisarenkoプロットを更新します。
1343        これにより、フィッティングプロセスのリアルタイムな監視が可能になります。
1344    引数:
1345        :param fig: matplotlibのFigureオブジェクト。
1346        :type fig: matplotlib.figure.Figure
1347        :param S2: 現在のS2値 (二乗和)。
1348        :type S2: float
1349    """
1350    global iter, graphupdateinterval
1351    global outputinterval
1352    global dos, mu
1353    global ysigma
1354    global method, tol, maxiter, h_diff
1355    global Nmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fit
1356
1357    if iter % outputinterval == 0:
1358        print("callback {:04d}: {:10.4g} {:10.4g}  S2={:10.6g}".format(iter, mobility.meff, mobility.l0, S2))
1359
1360    if iter % graphupdateinterval == 0:
1361        ndata = len(ysigma)
1362        ysigmafin = []
1363        ymufin    = []
1364        ySfin     = []
1365        for i in range(ndata):
1366            S      = yS[i]
1367            sigma  = ysigma[i]
1368            n      = yN[i]
1369            mu     = ymu[i]
1370
1371            if Nmin_fit is not None and n < Nmin_fit:
1372                continue
1373            if Nmax_fit is not None and Nmax_fit < n:
1374                continue
1375            if sigmamin_fit is not None and sigma < sigmamin_fit:
1376                continue
1377            if sigmamax_fit is not None and sigmamax_fit < sigma:
1378                continue
1379
1380            EFfin, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
1381            if not flag:
1382                app.terminate("Error in callback(): EF calculation did not converge. diffEF={}  eps={} in {} iteration"
1383                        .format(diffEF, epsEF, maxiter), usage = usage, pause = True)
1384
1385            sigmafin, nfin, mufin, tau_avg, Sfin, kappa, kappa_tot, L, PF, ZT, inf \
1386                = dos.cal_transport_properteis(T0, EFfin, mobility, klatt, validate_error_str = 'Error in fit()')
1387        
1388            ysigmafin.append(sigmafin)
1389            ymufin.append(mufin)
1390            ySfin.append(Sfin)
1391
1392        plot(fig, yN, ysigma, ymu, yS, ysigmaini, ymuini, ySini, ysigmafin, ymufin, ySfin)
1393
1394    iter += 1
1395
1396def callback_S(fig, xk):
1397    """
1398    概要:
1399        Seebeck係数のフィッティングのためのコールバック関数です。
1400    詳細説明:
1401        calS2_Sを呼び出してS2を計算し、現在のパラメータを更新してから、
1402        共通のコールバック関数callbackを呼び出して進行状況を表示します。
1403    引数:
1404        :param fig: matplotlibのFigureオブジェクト。
1405        :type fig: matplotlib.figure.Figure
1406        :param xk: 現在の最適化パラメータ (meff)。
1407        :type xk: numpy.ndarray
1408    """
1409    S2 = calS2_S(xk)
1410    ainew = [xk[0]]
1411    set_parameters(ainew)
1412
1413    callback(fig, S2)
1414
1415def callback_sigma(fig, xk):
1416    """
1417    概要:
1418        電気伝導率のフィッティングのためのコールバック関数です。
1419    詳細説明:
1420        calS2_sigmaを呼び出してS2を計算し、現在のパラメータを更新してから、
1421        共通のコールバック関数callbackを呼び出して進行状況を表示します。
1422    引数:
1423        :param fig: matplotlibのFigureオブジェクト。
1424        :type fig: matplotlib.figure.Figure
1425        :param xk: 現在の最適化パラメータ (l0)。
1426        :type xk: numpy.ndarray
1427    """
1428    global iter, graphupdateinterval
1429    global outputinterval
1430    global dos, mu
1431    global ysigma
1432    
1433    S2 = calS2_sigma(xk)
1434    ainew = [mobility.meff, xk[0]]
1435    set_parameters(ainew)
1436
1437    callback(fig, S2)
1438
1439def construct_lists(label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu):
1440    """
1441    概要:
1442        入力データのリストを正規化または補完します。
1443    詳細説明:
1444        sampleラベルまたはxsampleがNoneの場合、連番で補完します。
1445        sigma、N、muのいずれかがNoneの場合、他の利用可能なデータから計算を試みます。
1446    引数:
1447        :param label_sample: サンプルデータのラベル。
1448        :type label_sample: str or None
1449        :param xsample: サンプルデータのリスト。
1450        :type xsample: list
1451        :param label_S: Seebeck係数のラベル。
1452        :type label_S: str
1453        :param yS: Seebeck係数のリスト。
1454        :type yS: list
1455        :param label_sigma: 電気伝導率のラベル。
1456        :type label_sigma: str or None
1457        :param ysigma: 電気伝導率のリスト。
1458        :type ysigma: list
1459        :param label_N: キャリア濃度のラベル。
1460        :type label_N: str or None
1461        :param yN: キャリア濃度のリスト。
1462        :type yN: list
1463        :param label_mu: 移動度のラベル。
1464        :type label_mu: str or None
1465        :param ymu: 移動度のリスト。
1466        :type ymu: list
1467    戻り値:
1468        :returns: 正規化または補完されたラベルとデータリストのタプル。
1469        :rtype: tuple
1470    """
1471    ndata = len(yS)
1472    if label_sample is None:
1473        label_sample = 'sample'
1474        xsample = []
1475        for i in range(ndata):
1476            xsample.append(i + 1)
1477    if label_sigma is None:
1478        ysigma = []
1479        for i in range(ndata):
1480            ysigma.append(None)
1481    if label_N is None:
1482        yN = []
1483        for i in range(ndata):
1484            if ymu is not None:
1485                yN.append(None)
1486            else:
1487                yN.append(ysigma[i] / e / ymu[i])
1488    if label_mu is None:
1489        ymu = []
1490        for i in range(ndata):
1491            if yN[i] is not None:
1492                ymu.append(None)
1493            else:
1494                ymu.append(ysigma[i] / e / yN[i])
1495
1496    return label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu 
1497
1498def fit():
1499    """
1500    概要:
1501        実験データに対して熱電パラメータのフィッティングを行います。
1502    詳細説明:
1503        入力ファイルからSeebeck係数、電気伝導率、キャリア濃度、移動度のデータを読み込みます。
1504        指定された温度T0と初期パラメータ(meff, rfac, l0)に基づいて初期計算を行います。
1505        scipy.optimize.minimize を使用して、主に有効質量(meff)と平均自由行程プレファクタ(l0)をフィッティングし、
1506        計算値が実験データに最もよく合うように最適化します。
1507        フィッティング結果はパラメータファイルとExcelファイルに保存され、グラフで可視化されます。
1508        フィッティング範囲はNmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fitで指定できます。
1509    """
1510    global fig
1511    global label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu
1512    global ysigmaini, ymuini, ySini
1513    global method, tol, maxiter, h_diff
1514    global Nmin_fit, Nmax_fit, sigmamin_fit, sigmamax_fit
1515    
1516    print("infile     :", infile)
1517    print("  T_label    : ", T_label)
1518    print("  S_label    : ", S_label)
1519    print("  n_label    : ", n_label)
1520    print("  mu_label   : ", mu_label)
1521    print("  sigma_label: ", sigma_label)
1522    print("")
1523    print(f"T0: {T0} K")
1524    print("Carrier:")
1525    print("  q={} e".format(mobility.charge))
1526    print("  meff={}me".format(dos.meeff))
1527    dos.NC  = meff2NC_FEA(dos.meeff, T0)
1528    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
1529    print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
1530    print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
1531    print("  scattering factor (non-degenerated)     r={}".format(mobility.rfac))
1532    print("  mean free path prefactor: l0={} m".format(mobility.l0))
1533    print("")
1534    print( "Fitting condition")
1535    print(f"  method : {method}")
1536    print(f"  tol    : {tol}")
1537    print(f"  maxiter: {maxiter}")
1538    print(f"  h_diff: {h_diff}")
1539    print(f"  outputinterval: {outputinterval}")
1540    print(f"  N range    : {Nmin_fit} - {Nmax_fit} cm^-3")
1541    print(f"  sigma range: {sigmamin_fit} - {sigmamax_fit} S/cm")
1542    def set_default_range(var):
1543        if var == '' or var == '*':
1544            return None
1545        return pfloat(var, defval = None)
1546
1547    Nmin_fit = set_default_range(Nmin_fit)
1548    Nmax_fit = set_default_range(Nmax_fit)
1549    sigmamin_fit = set_default_range(sigmamin_fit)
1550    sigmamax_fit = set_default_range(sigmamax_fit)
1551
1552    if '***' in method:
1553        app.terminate(f"***Error: Choose method", pause = True)
1554
1555    if '***' in S_label:
1556        app.terminate(f"***Error: Choose S_label", pause = True)
1557    if '***' in n_label:
1558        app.terminate(f"***Error: Choose n_label", pause = True)
1559    if '***' in mu_label:
1560        app.terminate(f"***Error: Choose mu_label", pause = True)
1561
1562    if 'uV' in S_label or 'micro' in S_label or 'mV' in S_label:
1563        app.terminate(f"***Error: The unit of S must be 'V/K' but given by [{S_label}]", pause = True)
1564
1565    print("")
1566    print("Read S data from {}".format(infile))
1567
1568    datafile = tkVariousData(infile)
1569    labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = usage)
1570    label_S, yS           = datafile.FindDataArray(S_label, flag = 'i')
1571# Convert V/K to uV/K
1572    for i in range(len(yS)):
1573        yS[i] *= 1.0e6 
1574    label_sigma, ysigma   = datafile.FindDataArray(sigma_label, flag = 'i')
1575    label_N, yN           = datafile.FindDataArray(n_label, flag = 'i')
1576    label_mu, ymu         = datafile.FindDataArray(mu_label, flag = 'i')
1577    label_sample, xsample = labels[0], datalist[0]
1578
1579    ndata = len(yS)
1580    print("ndata(all): ", ndata)
1581    ndata_fit = 0
1582    print("data to be fitted")
1583    print(f"{'sample':15}  {'n(cm-3)':12}  {'mu(cm2/Vs)':10}  {'sigma(S/cm)':10}  {'S(uV/K)':10}")
1584    for i in range(ndata):
1585        sigma = ysigma[i]
1586        n = yN[i]
1587        if Nmin_fit is not None and n < Nmin_fit:
1588            continue
1589        if Nmax_fit is not None and Nmax_fit < n:
1590            continue
1591        if sigmamin_fit is not None and sigma < sigmamin_fit:
1592            continue
1593        if sigmamax_fit is not None and sigmamax_fit < sigma:
1594            continue
1595
1596        print(f"{xsample[i]:15}  {n:12.4g}  {ymu[i]:10.4g}  {sigma:10.4g}  {yS[i]:10.4g}")
1597
1598    print("")
1599    print("Calculate initial values")
1600    print_parameters()
1601    yEFini    = []
1602    ysigmaini = []
1603    ymuini    = []
1604    ySini     = []
1605    print("{:>8} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12} {:>12}"
1606            .format("EF(eV)", "Ne(cm-3)", "sigma(S/cm)", "mu(cm2/Vs)", "S(uV/K)", "Ne,ini", "sigma,ini", "mu,ini", "S,ini"))
1607    for i in range(ndata):
1608        sample = xsample[i]
1609        S      = yS[i]
1610        sigma  = ysigma[i]
1611        n      = yN[i]
1612        mu     = ymu[i]
1613
1614        EFini, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
1615        if not flag:
1616            app.terminate("Error in fit(): EF calculation did not converge. diffEF={}  eps={} in {} iteration".format(diffEF, epsEF, maxiter),
1617                        usage = usage, pause = True)
1618
1619#        x = (EF - doc.EC) * e / kB / T0
1620        sigmaini, nini, muini, tau_avg, Sini, kappa, kappa_tot, L, PF, ZT, inf = \
1621                dos.cal_transport_properteis(T0, EFini, mobility, klatt, validate_error_str = 'Error in fit()')
1622        
1623        yEFini.append(EFini)
1624        ysigmaini.append(sigmaini)
1625        ymuini.append(muini)
1626        ySini.append(Sini)
1627
1628        print("{:8.3g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g}"
1629                    .format(EFini, n, sigma, mu, S, nini, sigmaini, muini, Sini))
1630
1631#=============================
1632# グラフの表示
1633#=============================
1634    print("")
1635    print("plot")
1636    fig = plt.figure(figsize = figsize_sim)
1637
1638    plot(fig, yN, ysigma, ymu, yS, ysigmaini, ymuini, ySini)
1639
1640#=============================
1641# Optimization
1642#=============================
1643    print("")
1644    print("Nonlinear least-squares fitting for effective mass by ", method)
1645    print("  tol=", tol)
1646    ai = [mobility.meff]
1647    ret = minimize(calS2_S, ai, method = method, jac = diff1_S, tol = tol, callback = lambda xk: callback_S(fig, xk),
1648                options = {'maxiter':maxiter, "disp":True})
1649    if method == 'nelder-mead':
1650        simplex = ret['final_simplex']
1651        ai = simplex[0][0]    
1652        fmin = ret['fun']
1653    aifin = [ai[0], mobility.l0] #ai0[1]]
1654    set_parameters(aifin)
1655    print("Optimized at S2={:12.6g}:".format(fmin))
1656    print_parameters(parameter_list())
1657
1658    print("")
1659    print("Nonlinear least-squares fitting for l0 by ", method)
1660    print("  tol=", tol)
1661    ai_sigma = [mobility.l0]
1662    print("meff(ini)=", mobility.meff)
1663#    ret = minimize(calS2_S, ai_sigma, method = method, jac = diff1_sigma, tol = tol, callback = lambda xk: callback_sigma(fig, xk),
1664    ret = minimize(calS2_sigma, ai_sigma, method = method, jac = diff1_sigma, tol = tol, callback = lambda xk: callback_sigma(fig, xk),
1665                options = {'maxiter':maxiter, "disp":True})
1666    if method == 'nelder-mead':
1667        simplex = ret['final_simplex']
1668        ai_sigma = simplex[0][0]    
1669        fmin = ret['fun']
1670    aifin = [mobility.meff, ai_sigma[0]]
1671    set_parameters(aifin)
1672    print("Optimized at S2={:12.6g}:".format(fmin))
1673#    print_parameters()
1674    print_parameters(parameter_list())
1675    
1676    print("")
1677    print("Calculate final values")
1678    print_parameters(parameter_list())
1679    yEFfin    = []
1680    ysigmafin = []
1681    ySfin     = []
1682    print("{:>8} {:>12} {:>12} {:>12} {:>12} {:>12}"
1683            .format("EF(eV)", "Ne(cm-3)", "sigma(S/cm)", "mu(cm2/Vs)", "S(uV/K)", "S,fin"))
1684    for i in range(ndata):
1685        sample = xsample[i]
1686        S      = yS[i]
1687        sigma  = ysigma[i]
1688        n      = yN[i]
1689        mu     = ymu[i]
1690
1691        EFfin, diffEF, flag = dos.EF_from_electrondensity(n, T0, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
1692        if not flag:
1693            app.terminate("Error in fit(): EF calculation did not converge. diffEF={}  eps={} in {} iteration".format(diffEF, epsEF, maxiter),
1694                        usage = usage, pause = True)
1695
1696        sigmafin, nfin, mufin, tau_avg, Sfin, kappa, kappa_tot, L, PF, ZT, inf \
1697            = dos.cal_transport_properteis(T0, EFfin, mobility, klatt, validate_error_str = 'Error in fit()')
1698        
1699        mufin = sigmafin / sigma / mufin
1700        sigmafin = e * nfin * mufin
1701        
1702        yEFfin.append(EFfin)
1703        ysigmafin.append(sigmafin)
1704        ySfin.append(Sfin)
1705        print("{:8.3g} {:12.4g} {:12.4g} {:12.4g} {:12.4g} {:12.4g}"
1706                    .format(EFfin, n, sigma, mu, S, Sfin))
1707
1708    print("")
1709    print("Save final parameters to [{}]".format(parameterfile))
1710    save_parameterfile(ai = aifin, S2 = fmin)
1711
1712    print("Save fitting data to [{}]".format(outxlsfile))
1713    tkVariousData().to_excel(outxlsfile, ["EF(eV)", "Ne(cm-3)", "sigma(S/cm)", "mu(cm2/Vs)", "S(uV/K)", "S,fin"], 
1714            [yEFfin, yN, ysigma, ymu, yS, ySfin])
1715
1716    app.terminate("", pause = True)
1717
1718def sim():
1719    """
1720    概要:
1721        熱電特性とHall特性のシミュレーションを行い、結果をファイルとグラフに出力します。
1722    詳細説明:
1723        入力ファイルからSeebeck係数、電気伝導率、キャリア濃度、移動度、温度のデータを読み込みます。
1724        指定されたパラメータ (T0, meff, rfac, l0) を用いて、実験データ点における加重移動度を計算します。
1725        さらに、指定されたキャリア濃度範囲NminからNmaxまでを掃引し、
1726        Seebeck係数、電気伝導率、移動度、Hall係数、Hall因子、Hall移動度などの特性を計算します。
1727        計算結果はパラメータファイルとExcelファイルに保存され、PisarenkoプロットとJonkerプロットを含むグラフで可視化されます。
1728    """
1729    global label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu
1730    global Nmin, Nmax
1731    
1732    print("")
1733    print("infile     :", infile)
1734    print("  T_label    : ", T_label)
1735    print("  S_label    : ", S_label)
1736    print("  n_label    : ", n_label)
1737    print("  mu_label   : ", mu_label)
1738    print("  sigma_label: ", sigma_label)
1739    print(f"T0        : {T0} K")
1740    print("Carrier:")
1741    print("  q={} e".format(mobility.charge))
1742    print("  meff={}me".format(dos.meeff))
1743    dos.NC  = meff2NC_FEA(dos.meeff, T0)
1744    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
1745    if mobility.charge < 0.0:
1746        print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
1747        print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
1748    else:
1749        print("  NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
1750        print("  DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
1751    print("  scattering factor (non-degenerated)     r={}".format(mobility.rfac))
1752    print("  mean free path prefactor: l0={} m".format(mobility.l0))
1753
1754    if '***' in S_label:
1755        app.terminate(f"***Error: Choose S_label", pause = True)
1756    if '***' in n_label:
1757        app.terminate(f"***Error: Choose n_label", pause = True)
1758    if '***' in mu_label:
1759        app.terminate(f"***Error: Choose mu_label", pause = True)
1760
1761    if 'uV' in S_label or 'micro' in S_label or 'mV' in S_label:
1762        print("")
1763        print(f"***Error: The unit of S must be 'V/K' but given by [{S_label}]")
1764        app.terminate(pause = True)
1765
1766    print("")
1767    print("Read S data from {}".format(infile))
1768    datafile = tkVariousData(infile)
1769    labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = usage)
1770    label_S, yS           = datafile.FindDataArray(S_label, flag = 'i')
1771# Convert V/K to uV/K
1772    for i in range(len(yS)):
1773        yS[i] *= 1.0e6 
1774    label_sigma, ysigma   = datafile.FindDataArray(sigma_label, flag = 'i')
1775    label_N, yN           = datafile.FindDataArray(n_label, flag = 'i')
1776    label_mu, ymu         = datafile.FindDataArray(mu_label, flag = 'i')
1777    label_sample, xsample = labels[0], datalist[0]
1778    if T_label == '' or '***' in T_label:
1779        label_T = 'T(K)'
1780        yT = []
1781        for i in range(len(xsample)):
1782            yT.append(T0)
1783    else:
1784        label_T, yT           = datafile.FindDataArray(T_label, flag = 'i')
1785
1786#    print("x=", xsample)
1787    ndata = len(yS)
1788    print("ndata(all): ", ndata)
1789#    print("xsample=", label_sample, xsample)
1790#    print("yS=", label_S, yS)
1791#    print("ysigma=", label_sigma, ysigma)
1792    print(f"{'sample':15}  {'n(cm-3)':12}  {'mu(cm2/Vs)':10}  {'sigma(S/cm)':10}  {'S(uV/K)':10}")
1793    for i in range(ndata):
1794        sigma = ysigma[i]
1795        n = yN[i]
1796        print(f"{xsample[i]:15}  {n:12.4g}  {ymu[i]:10.4g}  {sigma:10.4g}  {yS[i]:10.4g}")
1797
1798
1799    print("")
1800    print("Weighted mobility")
1801    ysigma_w = []
1802    yn_w     = []
1803    ymu_w    = []
1804    ymu_w0   = []
1805    yEF      = []
1806    if yN[0] is not None:
1807        print("{:>10}\t{:>8}\t{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>15}\t{:>15}"
1808            .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)"))
1809    else:
1810        print("{:>10}\t{:>10}\t{:>10}\t{:>10}".format("sample", "EF(eV)", "S(uV/K)", "sigma(S/cm)"))
1811    for i in range(ndata):
1812        sample = xsample[i]
1813        S      = yS[i]
1814        sigma  = ysigma[i]
1815        n      = yN[i]
1816        mu     = ymu[i]
1817        T      = yT[i]
1818        print(f"sample={sample} S={S} uV/K sigma={sigma} S/cm")
1819        mu_w, n_w, sign, carriertype = weighted_mobility(sigma / 0.01, S * 1.0e-6, T0)
1820        mu_w   *= 1.0 * 1.0e4                       # cm2/Vs
1821        mu_w0  = mu_w * pow(mobility.meff, -1.5)
1822        n_w    *= 1.0e-6                            # cm^-3
1823        sigma_w = e * n_w * mu_w
1824
1825        if n is not None:
1826            EF, diffEF, flag = dos.EF_from_electrondensity(n, T, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
1827        else:
1828            EF, diffEF, flag = dos.EF_from_electronSeebeck(S, T, mobility, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100)
1829        if not flag:
1830            app.terminate("Error in sim(): EF calculation did not converge. diffEF={}  eps={} in {} iteration".format(diffEF, epsEF, maxiter),
1831                                usage = usage, pause = True)
1832#        print("n=", n, EF)       
1833
1834        ysigma_w.append(sigma_w)
1835        yn_w.append(n_w)
1836        ymu_w.append(mu_w)
1837        ymu_w0.append(mu_w0)
1838        yEF.append(EF)
1839        if n is not None:
1840            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}"
1841                .format(sample, T, EF, S, sigma, mu, n, sigma_w, mu_w, mu_w0, n_w))
1842        else:
1843            print("{:10.4g}\t{:8.3g}\t{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, T, EF, S, sigma))
1844
1845    Nmin_data = min(yN)
1846    Nmax_data = max(yN)
1847    if Nmin_data < Nmin:
1848        Nmin = Nmin_data
1849    if Nmax < Nmax_data:
1850        Nmax = Nmax_data
1851
1852    lnNmin = log(Nmin)
1853    lnNmax = log(Nmax)
1854    lnNstep = (lnNmax - lnNmin) / (nN - 1)
1855    print("")
1856    print("N range: {:10.4g} - {:10.4g} cm-3, {:10.4g} step in ln(N), nN={}".format(Nmin, Nmax, lnNstep, nN))
1857    xNsim        = []
1858    xsigmasim    = []
1859    yScal        = []
1860    ySsim_nondeg = []
1861    ySsim_deg    = []
1862    print("{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}"
1863            .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)"))
1864    for iN in range(nN):
1865        lnN = lnNmin + lnNstep * iN
1866        N   = exp(lnN)
1867#        print("N=", N, dos.NC)
1868        if N < dos.NC:
1869            EF = kB * T0 * log(N / dos.NC) / e
1870        else:
1871            EF = BMShift_FEA(N, dos.DC0)
1872
1873        sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
1874            = dos.cal_transport_properteis(T0, EF, mobility, klatt, validate_error_str = 'Error in sim()')
1875
1876        Sndeg = dos.cal_S_nondegenerated_from_Ne(n, mobility.rfac, charge = mobility.charge) * 1.0e6    # uV/K
1877        Sdeg  = dos.cal_S_degenerated_from_Ne(T0, n, mobility.rfac, charge = mobility.charge) * 1.0e6  # uV/K
1878
1879        xNsim.append(n)
1880        xsigmasim.append(sigma)
1881        yScal.append(S)
1882        ySsim_nondeg.append(Sndeg)
1883        ySsim_deg.append(Sdeg)
1884        print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}"
1885                .format(N, EF, S, Sndeg, Sdeg, sigma, n, mu))
1886
1887    print("")
1888    print("Save parameters to [{}]".format(parameterfile))
1889    save_parameterfile(S2 = '')
1890
1891    print("")
1892    print("Save data to [{}]".format(outxlsfile))
1893    if yN[0] is not None:
1894        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"],
1895                         [xsample, yS, ysigma, ymu, yN, ysigma_w, ymu_w, ymu_w0, yn_w])
1896    else:
1897        tkVariousData().to_excel(outxlsfile, ["sample", 'S(uV/K)', 'sigma(S/cm)'], [xsample, yS, ysigma])
1898
1899#=============================
1900# グラフの表示
1901#=============================
1902    print("")
1903    fig = plt.figure(figsize = figsize_sim)
1904
1905    ax1  = fig.add_subplot(2, 2, 1)
1906    ax1b = ax1.twinx()
1907    ax2  = fig.add_subplot(2, 2, 2)
1908    ax4  = fig.add_subplot(2, 2, 3)
1909    ax4b = ax4.twinx()
1910    ax3  = fig.add_subplot(2, 2, 4)
1911
1912    maxS = max(yS)
1913    minS = min(yS)
1914    if maxS > 0.0:
1915        maxS *= 2.0
1916    else:
1917        maxS = 0.0
1918    if minS > 0.0:
1919        minS = 0.0
1920    else:
1921        minS *= 2.0
1922
1923    ins1 = ax1.plot( xsample, yS,       label = 'S',          linestyle = '-',      linewidth = 1.0, color = 'red',   marker = 'o', markersize = 5.0)
1924    ins2 = ax1b.plot(xsample, ysigma,   label = r'$\sigma$',   linestyle = '-',      linewidth = 1.0, color = 'blue',  marker = '^', markersize = 10.0)
1925    ins3 = ax1b.plot(xsample, ysigma_w, label = r'$\sigma_w$', linestyle = 'dashed', linewidth = 0.5, color = 'green', marker = '+', markersize = 15.0, markerfacecolor = 'green')
1926    ax1.set_xlabel("sample", fontsize = fontsize)
1927    ax1.set_ylabel(r"S (K/$\mu$V)", fontsize = fontsize)
1928    ax1b.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
1929#    ax1.set_ylim([minS, maxS])
1930#    h1a, l1a = ax1.get_legend_handles_labels()
1931#    h1b, l1b = ax1b.get_legend_handles_labels()
1932#    ax1.legend(h1a + h1b, l1a + h1b, fontsize = legend_fontsize, loc = 0)
1933    ins = ins1 + ins2 + ins3
1934    ax1.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
1935#    ax1.legend(loc = 0)
1936#    ax1b.legend(loc = 0)
1937    ax1.tick_params(labelsize = fontsize)
1938    ax1b.tick_params(labelsize = fontsize)
1939
1940    if yN[0] is not None:
1941        ins1 = ax4.plot( xsample, ymu,    label = r'$\mu$',     linestyle = '-',      color = 'red',  linewidth = 0.5, marker = 'o', markersize = 10.0)
1942        ins2 = ax4.plot( xsample, ymu_w,  label = r'$\mu_w$',   linestyle = 'dashed', color = 'red',  linewidth = 0.5, marker = '^', markersize = 5.0)
1943        ins3 = ax4.plot( xsample, ymu_w0, label = r'$\mu_w^0$', linestyle = 'dashed', color = 'red',  linewidth = 0.5, marker = '+', markersize = 5.0)
1944        ins4 = ax4b.plot(xsample, yN,     label = '$N_e$',     linestyle = '-',      color = 'blue', linewidth = 0.5, marker = 'o', markersize = 10.0)
1945        ins5 = ax4b.plot(xsample, yn_w,   label = '$N_{e,w}$', linestyle = 'dashed', color = 'blue', linewidth = 0.5, marker = '^', markersize = 5.0)
1946        ax4.set_xlabel("sample", fontsize = fontsize)
1947        ax4.set_ylabel(r"$\mu$ (cm$^2$/V/s)", fontsize = fontsize)
1948        ax4b.set_ylabel(r"$N_e$ (cm$^{-3}$)", fontsize = fontsize)
1949#    ax4.set_ylim([minS, maxS])
1950        ax4b.set_yscale('log')
1951        ins = ins1 + ins2 + ins3 + ins4 + ins5  
1952        ax4.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
1953        ax4.tick_params(labelsize = fontsize)
1954        ax4b.tick_params(labelsize = fontsize)
1955
1956    ax2.plot(ysigma,    yS,           label = 'Jonker plot',  linestyle = '',       color = 'red',   marker = 'o', markersize = 5.0)
1957    ax2.plot(xsigmasim, yScal,        label = 'S,cal',        linestyle = '-',      color = 'red',   linewidth = 1.0)
1958    ax2.plot(xsigmasim, ySsim_nondeg, label = 'S,sim,nondeg', linestyle = 'dashed', color = 'blue',  linewidth = 0.5)
1959    ax2.plot(xsigmasim, ySsim_deg,    label = 'S,sim,deg',    linestyle = 'dashed', color = 'green', linewidth = 0.5)
1960    ax2.set_xlabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
1961    ax2.set_ylabel(r"S ($\mu$V/K)", fontsize = fontsize)
1962    ax2.set_xscale('log')
1963    ax2.set_ylim([minS, maxS])
1964    ax2.legend(fontsize = legend_fontsize, loc = 'best')
1965    ax2.tick_params(labelsize = fontsize)
1966
1967    if yN[0] is not None:
1968        ax3.plot(yN,    yS,           label = 'Pisarenko plot', linestyle = '', marker = 'o', markersize = 5.0)
1969    ax3.plot(xNsim, yScal,        label = 'S,cal',          linestyle = '-', color = 'red', linewidth = 1.0)
1970    ax3.plot(xNsim, ySsim_nondeg, label = 'S,sim,nondeg',   linestyle = 'dashed', color = 'blue', linewidth = 0.5)
1971    ax3.plot(xNsim, ySsim_deg,    label = 'S,sim,deg',      linestyle = 'dashed', color = 'green', linewidth = 0.5)
1972    ax3.set_xlabel(r"$N$ (cm$^{-3}$)", fontsize = fontsize)
1973    ax3.set_ylabel(r"S ($\mu$V/K)", fontsize = fontsize)
1974    ax3.set_xscale('log')
1975    ax3.set_ylim([minS, maxS])
1976    ax3.legend(fontsize = legend_fontsize, loc = 'best')
1977    ax3.tick_params(labelsize = fontsize)
1978
1979    plt.tight_layout()
1980    plt.pause(0.1)
1981
1982    app.terminate("", pause = True)
1983
1984def init():
1985    """
1986    概要:
1987        入力データファイルを読み込み、初期の熱電パラメータファイルを作成します。
1988    詳細説明:
1989        指定された入力ファイルからSeebeck係数、電気伝導率、キャリア濃度、移動度のデータを読み込み、
1990        基本的な情報を表示します。
1991        Seebeck係数の符号に基づいてキャリアの種類 (正孔または電子) を自動的に判断し、
1992        mobility.chargeを設定します。
1993        現在のパラメータはparameterfileに保存されます。
1994    """
1995    print("")
1996    print("Carrier:")
1997    print("  q={} e".format(mobility.charge))
1998    print("  meff={}me".format(dos.meeff))
1999    dos.NC  = meff2NC_FEA(dos.meeff, T0)
2000    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
2001    if mobility.charge < 0.0:
2002        print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
2003        print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
2004    else:
2005        print("  NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
2006        print("  DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
2007    print("  scattering factor (non-degenerated)     r={}".format(mobility.rfac))
2008    print("  mean free path prefactor: l0={} m".format(mobility.l0))
2009
2010    if 'uV' in S_label or 'micro' in S_label or 'mV' in S_label:
2011        print("")
2012        print(f"***Error: The unit of S must be 'V/K' but given by [{S_label}]")
2013        app.terminate(pause = True)
2014
2015    print("")
2016    print("Read S data from {}".format(infile))
2017    label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu \
2018            = read_datafile(infile, usage = usage)
2019# Convert V/K to uV/K
2020    for i in range(len(yS)):
2021        yS[i] *= 1.0e6 
2022
2023    ndata = len(xsample)
2024    print("ndata: ", ndata)
2025
2026    if yN[0] is not None:
2027        print("{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}"
2028            .format("sample", "S(uV/K)", "sigma(S/cm)", "mu(cm2/Vs)", "Ne(cm-3)"))
2029    else:
2030        print("{:>10}\t{:>10}\t{:>10}".format("sample", "S(uV/K)", "sigma(S/cm)"))
2031    for i in range(ndata):
2032        sample = xsample[i]
2033        S      = yS[i]
2034        sigma  = ysigma[i]
2035        n      = yN[i]
2036        mu     = ymu[i]
2037        if n is not None:
2038            print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, S, sigma, mu, n))
2039        else:
2040            print("{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, S, sigma))
2041
2042    if yS[0] > 0.0:
2043        mobility.charge = 1.0
2044    elif yS[0] < 0.0:
2045        mobility.charge = -1.0
2046    for i in range(1, ndata):
2047        if yS[i] * mobility.charge < 0.0:
2048            app.terminate("Error in init(): S data include both positive and negative values: S[{}}={}, S[{}]={}".format(0, yS[0], i, yS[i]),
2049                                usage = usage, pause = True)
2050
2051    print("")
2052    print("Save parameters to [{}]".format(parameterfile))
2053    save_parameterfile(S2 = '')
2054
2055def calL():
2056    """
2057    概要:
2058        温度、Hallキャリア濃度、電気伝導率のデータから電子のローレンツ数と熱伝導率を計算します。
2059    詳細説明:
2060        入力ファイルから温度、Hallキャリア濃度、電気伝導率のデータを読み込みます。
2061        指定された有効質量meff、散乱因子rfac、平均自由行程プレファクタl0に基づいて、
2062        各データ点におけるフェルミ準位EF、電子のローレンツ数L、および電子熱伝導率kappa_eを計算します。
2063        結果はコンソールとExcelファイルに出力され、グラフで可視化されます。
2064    """
2065    print("")
2066    print("Carrier:")
2067    print("  q={} e".format(mobility.charge))
2068    print("  meff={}me".format(dos.meeff))
2069    dos.NC  = meff2NC_FEA(dos.meeff, T0)
2070    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
2071    if mobility.charge < 0.0:
2072        print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
2073        print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
2074    else:
2075        print("  NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
2076        print("  DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
2077    print("  scattering factor (non-degenerated)     r={}".format(mobility.rfac))
2078    print("  mean free path prefactor: l0={} m".format(mobility.l0))
2079
2080    print("")
2081    print("Read T, Ne, and sigma data from {}".format(infile))
2082    datafile = tkVariousData(infile)
2083    labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = usage)
2084#    print("labels: ", labels)
2085#    ncol  = len(datalist)
2086#    print("ncol: ", ncol)
2087    label_T, xT           = datafile.FindDataArray(r'^T[\s|\(|\[|$]', flag = 'i')
2088    label_sigma, ysigma   = datafile.FindDataArray(r'^sigma[$|\S]?.*$', flag = 'i')
2089    label_N, yN           = datafile.FindDataArray(r'^N[$|\S]?.*$', flag = 'i')
2090    ndata = len(xT)
2091    print("ndata: ", ndata)
2092#    print("xT=", label_T, xT)
2093#    print("ysigma=", label_sigma, ysigma)
2094#    print("yN=", label_N, yN)
2095
2096    print("Calculat EF, L and kappa,e")
2097    yEF     = []
2098    yL      = []
2099    ykappa  = []
2100    print("kappa,e=sigma,obs*L,cal*T")
2101    print("{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}"
2102        .format("T(K)", "EF(eV)", "sigma(S/cm)", "Ne(cm^-3)", "L(Wohm/K^2)", "kappa,e"))
2103    for i in range(len(xT)):
2104        T = xT[i]
2105        n = yN[i]
2106        s = ysigma[i]
2107
2108        if n is not None:
2109            EF, diffEF, flag = dos.EF_from_electrondensity(n, T, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
2110        if not flag:
2111            app.terminate("Error in calL(): EF calculation did not converge. diffEF={}  eps={} in {} iteration".format(diffEF, epsEF, maxiter),
2112                            usage = usage, pause = True)
2113#        print("n=", n, EF)       
2114
2115        sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
2116            = dos.cal_transport_properteis(T, EF, mobility, klatt, validate_error_str = 'Error in properteis()')
2117
2118        kappa = L * s / 0.01 * T
2119
2120        yEF.append(EF)
2121        yL.append(L)
2122        ykappa.append(kappa)
2123
2124        print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}"
2125            .format(xT[i], EF, ysigma[i], yN[i], L, kappa))
2126
2127    header, ext     = os.path.splitext(infile)
2128    filebody        = os.path.basename(header)
2129    outxlsLfile     = f'{filebody}-L-me{dos.meeff}-r{mobility.rfac}.xlsx'
2130    print("")
2131    print("Save data to [{}]".format(outxlsLfile))
2132    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)"], 
2133                          [xT, yEF, ysigma, yN, yL, ykappa])
2134
2135#=============================
2136# グラフの表示
2137#=============================
2138    print("")
2139    fig = plt.figure(figsize = figsize_sim)
2140
2141    if min(xT) == max(xT):
2142        xX = yN
2143        xlabel = '$N$ (cm$^{-3}$)'
2144        xscale = 'log'
2145    else:
2146        xX = xT
2147        xlabel = '$T$ (K)'
2148        xscale = 'linear'
2149
2150    axEF    = fig.add_subplot(2, 2, 1)
2151    axNe    = fig.add_subplot(2, 2, 2)
2152    axsigma = axNe.twinx()
2153    axL     = fig.add_subplot(2, 2, 3)
2154    axkappa = axL.twinx()
2155    axkappa = fig.add_subplot(2, 2, 4)
2156
2157    axEF.plot(xX, yEF, label = 'EF(cal)', linestyle = '-', linewidth = 1.0,  marker = 'o', markersize = 5.0)
2158    axEF.set_xlabel(xlabel, fontsize = fontsize)
2159    axEF.set_ylabel("$E_F$ (eV)", fontsize = fontsize)
2160    axEF.set_xscale(xscale)
2161    axEF.legend(fontsize = legend_fontsize, loc = 'best')
2162    axEF.tick_params(labelsize = fontsize)
2163
2164    ins1  = axNe.plot   (xX, yN,     label = '$N_e$(obs)',           linestyle = '-', linewidth = 1.0, color = 'red',    marker = 'o', markersize = 5.0)
2165    ins2  = axsigma.plot(xX, ysigma, label = r'$\sigma$(obs)',        linestyle = 'dashed', linewidth = 1.0, color = 'blue', marker = 's', markersize = 5.0)
2166    axNe.set_xlabel(xlabel, fontsize = fontsize)
2167    axNe.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
2168    axNe.set_xscale(xscale)
2169    axNe.set_yscale('log')
2170    axsigma.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
2171    axsigma.set_yscale('log')
2172    ins = ins1 + ins2
2173    axNe.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
2174    axNe.tick_params(labelsize = fontsize)
2175    axsigma.tick_params(labelsize = fontsize)
2176
2177    axL.plot(xX, yL, label = 'L(cal)', linestyle = '-', linewidth = 1.0,  marker = 'o', markersize = 5.0)
2178    axL.set_xlabel(xlabel, fontsize = fontsize)
2179    axL.set_ylabel(r"$L$ (W$\Omega$/K$^2$)", fontsize = fontsize)
2180    axL.set_xscale(xscale)
2181    axL.legend(fontsize = legend_fontsize, loc = 'best')
2182    axL.tick_params(labelsize = fontsize)
2183    
2184    axkappa.plot(xX, ykappa, label = r'$\kappa_e$(cal)', linestyle = '-', linewidth = 1.0,  marker = 'o', markersize = 5.0)
2185    axkappa.set_xlabel(xlabel, fontsize = fontsize)
2186    axkappa.set_ylabel(r"$\kappa_e$=$\sigma$$L$$T$ (W/m/K)", fontsize = fontsize)
2187    axkappa.set_xscale(xscale)
2188    axkappa.set_yscale('log')
2189    axkappa.legend(fontsize = legend_fontsize, loc = 'best')
2190    axkappa.tick_params(labelsize = fontsize)
2191
2192    plt.tight_layout()
2193    plt.pause(0.1)
2194
2195    app.terminate("", pause = True)
2196
2197def calF():
2198    """
2199    概要:
2200        温度、Hallキャリア濃度、Hall移動度のデータからHall因子、ドリフトキャリア濃度、ドリフト移動度を計算します。
2201    詳細説明:
2202        入力ファイルから温度、Hallキャリア濃度、Hall移動度、電気伝導率のデータを読み込みます。
2203        指定された有効質量meff、散乱因子rfac、平均自由行程プレファクタl0に基づいて、
2204        各データ点におけるフェルミ準位EF、Hall因子FH、ドリフトキャリア濃度ndrift、およびドリフト移動度mudriftを計算します。
2205        結果はコンソールとExcelファイルに出力され、グラフで可視化されます。
2206    """
2207    print("")
2208    print("infile     :", infile)
2209    print("T_label    : ", T_label)
2210    print("n_label    : ", n_label)
2211    print("mu_label   : ", mu_label)
2212    print("sigma_label: ", sigma_label)
2213    print("Carrier:")
2214    print("  q={} e".format(mobility.charge))
2215    print("  meff={}me".format(dos.meeff))
2216    dos.NC  = meff2NC_FEA(dos.meeff, T0)
2217    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
2218    if mobility.charge < 0.0:
2219        print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
2220        print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
2221    else:
2222        print("  NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
2223        print("  DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
2224    print("  scattering factor (non-degenerated)     r={}".format(mobility.rfac))
2225    print("  mean free path prefactor: l0={} m".format(mobility.l0))
2226
2227    print("")
2228    print("Read T, Ne, and sigma data from {}".format(infile))
2229    datafile = tkVariousData(infile)
2230    labels, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = usage)
2231    label_T, xT           = datafile.FindDataArray(T_label, flag = 'i')
2232    label_sigma, ysigma   = datafile.FindDataArray(sigma_label, flag = 'i')
2233    label_N, yN           = datafile.FindDataArray(n_label, flag = 'i')
2234    label_mu, ymu         = datafile.FindDataArray(mu_label, flag = 'i')
2235
2236    ndata = len(xT)
2237    print("ndata: ", ndata)
2238
2239    print("")
2240    print("Calculat EF, FH, mu_drift, and Ne")
2241    yEF      = []
2242    yFH      = []
2243    yRH      = []
2244    yndrift  = []
2245    ymudrift = []
2246    print("kappa,e=sigma,obs*L,cal*T")
2247    print("{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}"
2248        .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)"))
2249    for i in range(len(xT)):
2250        T = xT[i]
2251        n = yN[i]
2252        s = ysigma[i]
2253
2254        if n is not None:
2255            EF, diffEF, flag = dos.EF_from_electrondensity(n, T, EF0 = 0.0, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
2256        if not flag:
2257            app.terminate("Error in calF(): EF calculation did not converge. diffEF={}  eps={} in {} iteration".format(diffEF, epsEF, maxiter), 
2258                        usage = usage, pause = True)
2259
2260#        sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
2261#            = dos.cal_transport_properteis(T, EF, mobility, klatt, validate_error_str = 'Error in properteis()')
2262        inf = dos.cal_Hall_properteis(T0, EF, mobility, validate_error_str = 'Error in properteis()')
2263        FH       = inf["FH"]
2264        RH       = inf["RH"]
2265        RH0      = inf["RH0"]
2266        n_drift  = inf["nHall"]
2267        mu_drift = inf["muHall"]
2268
2269        yEF.append(EF)
2270        yFH.append(FH)
2271        yRH.append(RH)
2272        yndrift.append(yN[i] * FH)
2273        ymudrift.append(ymu[i] / FH)
2274
2275        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}"
2276            .format(xT[i], EF, ysigma[i], yN[i], ymu[i], yndrift[i], ymudrift[i], FH, RH))
2277
2278    header, ext     = os.path.splitext(infile)
2279    filebody        = os.path.basename(header)
2280    outxlsLfile     = f'{filebody}-FH-me{dos.meeff}-r{mobility.rfac}.xlsx'
2281    print("")
2282    print("Save data to [{}]".format(outxlsLfile))
2283    tkVariousData().to_excel(outxlsLfile, ["T(K)", "EF,cal(eV)", "sigma,obs(S/cm)", "Ne,Hall(obs)(cm-3)", "mu,Hall(obs)(cm2/Vs)", 
2284                           "Ne=Ne,Hall*FH(cm-3)", "mu,drift=mu,Hall/FH(cm2/Vs)", "FH", "RH(cm3/C)"], 
2285                          [xT, yEF, ysigma, yN, ymu, yndrift, ymudrift, yFH, yRH])
2286
2287#=============================
2288# グラフの表示
2289#=============================
2290    print("")
2291    fig = plt.figure(figsize = figsize_sim)
2292
2293    if min(xT) == max(xT):
2294        xX = yN
2295        xlabel = '$N$ (cm$^{-3}$)'
2296        xscale = 'log'
2297    else:
2298        xX = xT
2299        xlabel = '$T$ (K)'
2300        xscale = 'linear'
2301
2302    axEF    = fig.add_subplot(2, 3, 1)
2303    axsigma = fig.add_subplot(2, 3, 2)
2304    axNe    = fig.add_subplot(2, 3, 3)
2305    axmu    = fig.add_subplot(2, 3, 4)
2306    axFH    = fig.add_subplot(2, 3, 5)
2307    axRH    = fig.add_subplot(2, 3, 6)
2308
2309    axEF.plot(xX, yEF, label = 'EF(cal)', linestyle = '-', linewidth = 1.0,  marker = 'o', markersize = 5.0)
2310    axEF.set_xlabel(xlabel, fontsize = fontsize)
2311    axEF.set_ylabel("$E_F$ (eV)", fontsize = fontsize)
2312    axEF.set_xscale(xscale)
2313    axEF.legend(fontsize = legend_fontsize, loc = 'best')
2314    axEF.tick_params(labelsize = fontsize)
2315
2316    axsigma.plot(xX, ysigma, label = r'$\sigma$(obs)',        linestyle = 'dashed', linewidth = 1.0, color = 'blue', marker = 's', markersize = 5.0)
2317    axsigma.set_xlabel(xlabel, fontsize = fontsize)
2318    axsigma.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
2319    axsigma.set_xscale(xscale)
2320    axsigma.set_yscale('log')
2321    axsigma.legend(fontsize = legend_fontsize, loc = 'best')
2322    axsigma.tick_params(labelsize = fontsize)
2323
2324    axNe.plot   (xX, yN,      label = '$N_e$$_{,Hall}(obs)',   linestyle = '-', linewidth = 1.0, color = 'red',   marker = 'o', markersize = 5.0)
2325    axNe.plot   (xX, yndrift, label = '$N_e$$_{,drift}$(obs)', linestyle = '-', linewidth = 0.5, color = 'green', marker = '^', markersize = 3.0)
2326    axNe.set_xlabel(xlabel, fontsize = fontsize)
2327    axNe.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
2328    axNe.set_xscale(xscale)
2329    axNe.set_yscale('log')
2330    axNe.legend(fontsize = legend_fontsize, loc = 'best')
2331    axNe.tick_params(labelsize = fontsize)
2332
2333    axmu.plot(xX, ymu,      label = r'$\mu_{Hall}$(obs)',  linestyle = '-', linewidth = 1.0,  color = 'red',   marker = 'o', markersize = 5.0)
2334    axmu.plot(xX, ymudrift, label = r'$\mu_{drift}$(obs)', linestyle = '-', linewidth = 0.5,  color = 'green', marker = '^', markersize = 3.0)
2335    axmu.set_xlabel(xlabel, fontsize = fontsize)
2336    axmu.set_ylabel(r"$\mu$ (cm$^2$/Vs)", fontsize = fontsize)
2337    axmu.set_xscale(xscale)
2338    axmu.legend(fontsize = legend_fontsize, loc = 'best')
2339    axmu.tick_params(labelsize = fontsize)
2340
2341    axFH.plot(xX, yFH, label = '$F_{Hall}$', linestyle = '-', linewidth = 1.0,  color = 'blue', marker = 'o', markerfacecolor = 'blue', markersize = 3.0)
2342    axFH.set_xlabel(xlabel, fontsize = fontsize)
2343    axFH.set_ylabel("$F_{Hall}$", fontsize = fontsize)
2344    axFH.set_xscale(xscale)
2345    axFH.legend(fontsize = legend_fontsize, loc = 'best')
2346    axFH.tick_params(labelsize = fontsize)
2347
2348    axRH.plot(xX, yRH, label = '$R_H$',      linestyle = '-', linewidth = 1.0,  color = 'blue',  marker = 'o', markerfacecolor = 'blue', markersize = 3.0)
2349    axRH.set_xlabel(xlabel, fontsize = fontsize)
2350    axRH.set_ylabel("$R_H$ (cm$^3$/C)", fontsize = fontsize)
2351    axRH.set_xscale(xscale)
2352    axRH.set_yscale('log')
2353    axRH.legend(fontsize = legend_fontsize, loc = 'best')
2354    axRH.tick_params(labelsize = fontsize)
2355
2356    plt.tight_layout()
2357    plt.pause(0.1)
2358
2359    app.terminate("", pause = True)
2360
2361def T():
2362    """
2363    概要:
2364        温度とキャリア濃度に対する熱電特性の依存性をシミュレートします。
2365    詳細説明:
2366        入力ファイルからSeebeck係数、電気伝導率、キャリア濃度、移動度などのデータを読み込みます。
2367        指定された温度範囲 (TminからTmax) とキャリア濃度範囲 (NminからNmax) に対して、
2368        Seebeck係数、電気伝導率、移動度、電子熱伝導率、力率、ZT値、フェルミ準位を計算します。
2369        計算結果はコンソールと複数のExcelファイルに保存され、グラフで可視化されます。
2370    """
2371    global label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu
2372
2373    print("")
2374    print("Carrier:")
2375    print("  q={} e".format(mobility.charge))
2376    print("  meff={}me".format(dos.meeff))
2377    dos.NC  = meff2NC_FEA(dos.meeff, T0)
2378    dos.DC0 = meff2DC0_FEA(dos.meeff, T0)
2379    if mobility.charge < 0.0:
2380        print("  NC={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
2381        print("  DC={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
2382    else:
2383        print("  NV={:10.4g} cm^-3 at {} K".format(dos.NC, T0))
2384        print("  DV={:10.4g} cm^-3 at {} K".format(dos.DC0, T0))
2385    print("  scattering factor (non-degenerated)     r={}".format(mobility.rfac))
2386    print("  mean free path prefactor: l0={} m".format(mobility.l0))
2387
2388    print("")
2389    print("Read S data from {}".format(infile))
2390    label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu \
2391            = read_datafile(infile, usage = usage)
2392    label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu \
2393            = construct_lists(label_sample, xsample, label_S, yS, label_sigma, ysigma, label_N, yN, label_mu, ymu)
2394    print("x=", xsample)
2395    ndata = len(xsample)
2396    print("ndata: ", ndata)
2397
2398    if yN[0] is not None:
2399        print("{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}\t{:>10}"
2400            .format("sample", "EF(eV)", "S(uV/K)", "sigma(S/cm)", "mu(cm2/Vs)", "Ne(cm-3)"))
2401    else:
2402        print("{:>10}\t{:>10}\t{:>10}\t{:>10}".format("sample", "EF(eV)", "S(uV/K)", "sigma(S/cm)"))
2403    for i in range(ndata):
2404        sample = xsample[i]
2405        S      = yS[i]
2406        sigma  = ysigma[i]
2407        n      = yN[i]
2408        mu     = ymu[i]
2409        if n is not None:
2410            print("{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, S, sigma, mu, n))
2411        else:
2412            print("{:10.4g}\t{:10.4g}\t{:10.4g}".format(sample, S, sigma))
2413
2414    Nmin_data = min(yN)
2415    Nmax_data = max(yN)
2416    if Nmin_data < Nmin:
2417        Nmin = Nmin_data
2418    if Nmax < Nmax_data:
2419        Nmax = Nmax_data
2420    print("")
2421    print(f"Plot N range: {Nmin:12.4g} - {Nmax:12.4g}")
2422
2423    lnNmin = log(Nmin)
2424    lnNmax = log(Nmax)
2425    lnNstep = (lnNmax - lnNmin) / (nN - 1)
2426    print("")
2427    print("N range: {:10.4g} - {:10.4g} cm-3, {:10.4g} step in ln(N), nN={}".format(Nmin, Nmax, lnNstep, nN))
2428    xEFcal       = make_matrix2(nT, nN, 0.0)
2429    xNsim        = make_matrix2(nT, nN, 0.0)
2430    ysigmasim    = make_matrix2(nT, nN, 0.0)
2431    ymusim       = make_matrix2(nT, nN, 0.0)
2432    yScal        = make_matrix2(nT, nN, 0.0)
2433    ykappacal    = make_matrix2(nT, nN, 0.0)
2434    yPFcal       = make_matrix2(nT, nN, 0.0)
2435    yZTcal       = make_matrix2(nT, nN, 0.0)
2436    Tstep = (Tmax - Tmin) / (nT - 1)
2437    for iT in range(nT):
2438        T = Tmin + iT * Tstep
2439        dos.NC  = meff2NC_FEA(dos.meeff, T)
2440        dos.DC0 = meff2DC0_FEA(dos.meeff, T)
2441
2442        print("{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}"
2443                .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"))
2444        for iN in range(nN):
2445            lnN = lnNmin + lnNstep * iN
2446            N   = exp(lnN)
2447#           print("N=", N, dos.NC)
2448            if N < dos.NC:
2449                EF = kB * T0 * log(N / dos.NC) / e
2450            else:
2451                EF = BMShift_FEA(N, dos.DC0)
2452
2453            EF, diffEF, flag = dos.EF_from_electrondensity(N, T, EF0 = EF, dEF = 0.01, dump = 0.0, epsEF = epsEF, maxiter = 100) 
2454            if not flag:
2455               app.terminate("Error in sim(): EF calculation did not converge. diffEF={}  eps={} in {} iteration".format(diffEF, epsEF, maxiter),
2456                                usage = usage, pause = True)
2457
2458            sigma, n, mu, tau_avg, S, kappa, kappa_tot, L, PF, ZT, inf \
2459                = dos.cal_transport_properteis(T, EF, mobility, klatt, validate_error_str = 'Error in sim()')
2460
2461            xEFcal[iT][iN]     = EF
2462            xNsim[iT][iN]      = n
2463            ysigmasim[iT][iN]  = sigma
2464            ymusim[iT][iN]     = mu
2465            yScal[iT][iN]      = S
2466            ykappacal[iT][iN]  = kappa
2467            yPFcal[iT][iN]     = PF
2468            yZTcal[iT][iN]     = ZT
2469            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}"
2470                    .format(T, N, EF, S, sigma, n, mu, kappa, PF, ZT))
2471
2472    print("")
2473    print("Save parameters to [{}]".format(parameterfile))
2474    save_parameterfile(S2 = '')
2475
2476
2477    print("")
2478    labels = ["N(cm-3)"] + ["{:.4g} K".format(Tmin + iT * Tstep) for iT in range(nT)]
2479    header, ext     = os.path.splitext(infile)
2480    filebody        = os.path.basename(header)
2481    xlsfile      = filebody + '-EF.xlsx'.format(T)
2482    print("Save data to [{}]".format(xlsfile))
2483#    print("labels=", labels)
2484#    print("datalist=", [xNsim[0]])# + [xEFcal[iT] for iT in range(nT)])
2485    tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + xEFcal) #[xEFcal[iT] for iT in range(nT)])
2486    xlsfile   = filebody + '-sigma.xlsx'.format(T)
2487    print("Save data to [{}]".format(xlsfile))
2488    tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + ysigmasim)
2489    xlsfile   = filebody + '-mu.xlsx'.format(T)
2490    print("Save data to [{}]".format(xlsfile))
2491    tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + ymusim)
2492    xlsfile   = filebody + '-S.xlsx'.format(T)
2493    print("Save data to [{}]".format(xlsfile))
2494    tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + yScal)
2495    xlsfile   = filebody + '-kappa.xlsx'.format(T)
2496    print("Save data to [{}]".format(xlsfile))
2497    tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + ykappacal)
2498    xlsfile   = filebody + '-PF.xlsx'.format(T)
2499    print("Save data to [{}]".format(xlsfile))
2500    tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + yPFcal)
2501    xlsfile   = filebody + '-ZT.xlsx'.format(T)
2502    print("Save data to [{}]".format(xlsfile))
2503    tkVariousData().to_excel(xlsfile, labels, [xNsim[0]] + yZTcal)
2504
2505#=============================
2506# グラフの表示
2507#=============================
2508    print("")
2509    fig = plt.figure(figsize = figsize)
2510
2511    axSsobs    = fig.add_subplot(3, 3, 1)
2512    axSsobsb   = axSsobs.twinx()
2513    axmuNobs   = fig.add_subplot(3, 3, 2)
2514    axmuNobsb  = axmuNobs.twinx()
2515    axNs       = fig.add_subplot(3, 3, 3)
2516    axNmu      = fig.add_subplot(3, 3, 4)
2517    axNS       = fig.add_subplot(3, 3, 5)
2518    axNkappa   = fig.add_subplot(3, 3, 6)
2519    axNPF      = fig.add_subplot(3, 3, 7)
2520    axNZT      = fig.add_subplot(3, 3, 8)
2521    axNEF      = fig.add_subplot(3, 3, 9)
2522
2523    maxS = max(yS)
2524    minS = min(yS)
2525    if maxS > 0.0:
2526        maxS *= 1.2
2527    else:
2528        maxS = 0.0
2529    if minS > 0.0:
2530        minS = 0.0
2531    else:
2532        minS *= 1.2
2533
2534    ins1 = axSsobs.plot (xsample, yS,       label = 'S',          linestyle = '-',      linewidth = 1.0, color = 'red',   marker = 'o', markersize = 5.0)
2535    ins2 = axSsobsb.plot(xsample, ysigma,   label = r'$\sigma$',   linestyle = '-',      linewidth = 1.0, color = 'blue',  marker = '^', markersize = 10.0)
2536    axSsobs.set_xlabel("sample", fontsize = fontsize)
2537    axSsobs.set_ylabel(r"S (K/$\mu$V)", fontsize = fontsize)
2538    axSsobsb.set_ylabel(r"$\sigma$ (S/cm)", fontsize = fontsize)
2539    ins = ins1 + ins2 
2540    axSsobs.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'upper center') #loc = 'best')
2541    axSsobs.tick_params(labelsize = fontsize)
2542    axSsobsb.tick_params(labelsize = fontsize)
2543
2544    if yN[0] is not None:
2545        ins1 = axmuNobs.plot( xsample, ymu,    label = r'$\mu$',     linestyle = '-',      color = 'red',  linewidth = 0.5, marker = 'o', markersize = 10.0)
2546        ins2 = axmuNobsb.plot(xsample, yN,     label = '$N_e$',     linestyle = '-',      color = 'blue', linewidth = 0.5, marker = 'o', markersize = 10.0)
2547        axmuNobs.set_xlabel("sample", fontsize = fontsize)
2548        axmuNobs.set_ylabel(r"$\mu$ (cm$^2$/V/s)", fontsize = fontsize)
2549        axmuNobsb.set_ylabel("$N_e$ (cm$^{-3}$)", fontsize = fontsize)
2550        axmuNobsb.set_yscale('log')
2551        ins = ins1 + ins2
2552        axmuNobs.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
2553        axmuNobs.tick_params(labelsize = fontsize)
2554        axmuNobsb.tick_params(labelsize = fontsize)
2555
2556    if yN[0] is not None:
2557        axNs.plot( yN, ysigma, label = r'$\sigma_{,obs}$', linestyle = '', marker = 'o', markersize = 5.0)
2558        axNmu.plot(yN, ymu,    label = '$N_e$$_{,obs}$', linestyle = '', marker = 'o', markersize = 5.0)
2559        axNS.plot( yN, yS,     label = '$S_{obs}$', linestyle = '', marker = 'o', markersize = 5.0)
2560    for iT in range(nT):
2561        T= Tmin + iT * Tstep
2562        axNs.plot(    xNsim[iT], ysigmasim[iT], label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
2563        axNmu.plot(   xNsim[iT], ymusim[iT],    label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
2564        axNS.plot(    xNsim[iT], yScal[iT],     label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
2565        axNkappa.plot(xNsim[iT], ykappacal[iT], label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
2566        axNPF.plot(   xNsim[iT], yPFcal[iT],    label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
2567        axNZT.plot(   xNsim[iT], yZTcal[iT],    label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
2568        axNEF.plot(   xNsim[iT], xEFcal[iT],    label = '{:8.2f} K'.format(T), linestyle = '-', linewidth = 0.5)
2569    axNs.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
2570    axNmu.set_xlabel(r"$N$ (cm$^{-3}$)", fontsize = fontsize)
2571    axNS.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
2572    axNkappa.set_xlabel( "$N$ (cm$^{-3}$)", fontsize = fontsize)
2573    axNPF.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
2574    axNZT.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
2575    axNEF.set_xlabel( r"$N$ (cm$^{-3}$)", fontsize = fontsize)
2576    axNs.set_ylabel( r"$\sigma$ (S/cm)", fontsize = fontsize)
2577    axNmu.set_ylabel(r"$\mu$ (cm$^2$/V/s)", fontsize = fontsize)
2578    axNS.set_ylabel( "S (V/K)", fontsize = fontsize)
2579    axNkappa.set_ylabel( r"$\kappa$$_{,e}$ (W/m/K)", fontsize = fontsize)
2580    axNPF.set_ylabel( r"PF ($\mu$W/cm/K$^2$)", fontsize = fontsize)
2581    axNZT.set_ylabel( "ZT", fontsize = fontsize)
2582    axNEF.set_ylabel( "$E_F$", fontsize = fontsize)
2583    axNs.set_xscale('log')
2584    axNmu.set_xscale('log')
2585    axNS.set_xscale('log')
2586    axNkappa.set_xscale('log')
2587    axNPF.set_xscale('log')
2588    axNZT.set_xscale('log')
2589    axNEF.set_xscale('log')
2590    axNs.set_yscale('log')
2591    axNkappa.set_yscale('log')
2592    axNs.legend(fontsize = legend_fontsize, loc = 'best')
2593    axNmu.legend(fontsize = legend_fontsize, loc = 'best')
2594    axNS.legend(fontsize = legend_fontsize, loc = 'best')
2595    axNkappa.legend(fontsize = legend_fontsize, loc = 'best')
2596    axNPF.legend(fontsize = legend_fontsize, loc = 'best')
2597    axNZT.legend(fontsize = legend_fontsize, loc = 'best')
2598    axNEF.legend(fontsize = legend_fontsize, loc = 'best')
2599    
2600    plt.tight_layout()
2601    plt.pause(0.1)
2602
2603    app.terminate("", pause = True)
2604
2605
2606def main():
2607    """
2608    概要:
2609        スクリプトの主要な実行フローを制御します。
2610    詳細説明:
2611        アプリケーションの初期化、ログファイルのリダイレクト、コマンドライン引数に基づく変数の更新、
2612        パラメータファイルの読み書き、および選択されたモードに応じた各種計算関数の呼び出しを行います。
2613        電子構造および輸送パラメータの情報を出力します。
2614    """
2615    global app
2616    global parameterfile, parameterbkfile, outxlsfile, outxlsPropfile
2617    global outxlsfFDfile, outxlsFjfile
2618    
2619    app = tkApplication()
2620    logfile = app.replace_path(infile)
2621    print(f"Open logfile [{logfile}]")
2622    app.redirect(targets = ["stdout", logfile], mode = 'w')
2623
2624    updatevars()
2625
2626    parameterfile   = app.replace_path(infile, template = ["{dirname}", "{filebody}.in"])
2627    parameterbkfile = app.replace_path(infile, template = ["{dirname}", "{filebody}.in.bak"])
2628    outxlsfile      = app.replace_path(infile, template = ["{dirname}", "{filebody}-Seebeck-{mode}.xlsx"], ext_dict = {"mode": mode})
2629    outxlsPropfile  = app.replace_path(infile, template = ["{dirname}", "{filebody}-properties-{mode}.xlsx"], ext_dict = {"mode": mode})
2630
2631    outxlsfFDfile   = app.replace_path(infile, template = ["{dirname}", "fFD.xlsx"])
2632    outxlsFjfile    = app.replace_path(infile, template = ["{dirname}", "Fj.xlsx"])
2633
2634    print("")
2635    print("mode: ", mode)
2636    print("infile               : {}".format(infile))
2637    print("out xlsx file        : {}".format(outxlsfile))
2638    print("parameter file       : {}".format(parameterfile))
2639    print("parameter backup file: {}".format(parameterbkfile))
2640    print("")
2641
2642    print("")
2643    print("Read [{}]".format(parameterfile))
2644    read_parameters(parameterfile)
2645# check args again to use given parameters
2646    updatevars()
2647
2648    print("")
2649    print("Electronic structure")
2650#    print(f"Eg={dos.EC - dos.EV} eV")
2651#    print(f"EV={dos.EV} eV")
2652    print(f"EC={dos.EC} eV")
2653#    print(f"EA={dos.EA} eV")
2654#    print(f"NA={dos.NA} cm^-3")
2655#    print(f"ED={dos.ED} eV")
2656#    print(f"ND={dos.ND} cm^-3")
2657#    print(f"mh_eff={dos.mheff} me_0")
2658    print(f"me_eff={dos.meeff} me_0")
2659    print(f"Initial parameter")
2660    print(f"  EF0={dos.EF0} eV")
2661    print("Transport parameters")
2662    print(f"charge={mobility.charge} e")
2663    print(f"meff={mobility.meff}  me_0")
2664    print(f"rfac={mobility.rfac}")
2665    print(f"l0  ={mobility.l0} m")
2666
2667    if mode == 'basic':
2668        basic()
2669    elif mode == 'prop':
2670        properties()
2671    elif mode == 'Hall':
2672        Hall()
2673    elif mode == 'init':
2674        init()
2675    elif mode == 'sim':
2676        sim()
2677    elif mode == 'fit':
2678        fit()
2679    elif mode == 'T':
2680        T()
2681    elif mode == 'calL':
2682        calL()
2683    elif mode == 'calF':
2684        calF()
2685    else:
2686        app.terminate("Error in main: Invalide mode [{}]".format(mode), usage = usage, pause = True)
2687
2688if __name__ == "__main__":
2689    main()