import os
import sys
from math import exp, sqrt, log
import numpy as np
#from numpy import arange
from scipy import optimize # newton関数はscipy.optimizeモジュールに入っている
from scipy.interpolate import interp1d
from matplotlib import pyplot as plt
from tklib.tkfile import tkFile
from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg, save_csv
from tklib.tksci.tksci import pi, h, hbar,c, e, kB, me
from tklib.tktransport.tkTransport import integrate_Simpson_list, fe, fh, NC2meff_FEA, meff2DC0_FEA
from tklib.tktransport.tkDOS_VASP import integrate_Simpson, integrate_Simpson_list, tkDOS
"""
VASPのDOSCARに基づいてキャリア密度とフェルミ準位を計算するモジュール。
詳細説明:
VASPのファイルからDOSを読み込み、指定されたモード（me, T, EF）に応じて状態密度、
キャリア密度の温度依存性やフェルミ準位依存性を計算・プロットする。
関連リンク:
- :doc:`EF-T-semi_VASP_usage`
"""
#================================
# Global variables
#================================
Debug = 0
#mode: 'me', 'EF', 'T'
mode = 'me'
# For mode = 'me'. E range for fitting to me*(DOS)
dE_me_lsq = 0.5 # measured from EC/EV, eV
# EF range for mode = 'EF'
T0 = 300 # K
dEFmin = -0.5 # measured from EV, eV
dEFmax = 0.5 # measured from EC, eV
nEF = 50
EFmin = None
EFmax = None
EFstep = None
# Temperature range for mode = 'T'
Tmin = 300 # K
Tmax = 600
nT = 11
Tstep = None
# read from DOSCAR
dos = tkDOS()
file = "DOSCAR"
outcsvfile = None
# Flag to specify the data to determine band edges
#dos.data_for_bandedges = 'EIGENVAL'
dos.data_for_bandedges = 'DOSCAR'
dos.dEA = 0.05 # Acceptor level measured from EV, eV
dos.NA = 0.8e17
dos.dED = -0.28 # Donor level measured from EC, eV
dos.ND = 0.0e17
dos.EA = None # Acceptor states are not considered if EA = None
dos.ED = None # Donor states are not considered if ED = None
# other parameters
dos.EF0 = 0.0 # EF at charge neutral
dos.Egth = 1.0e18 # Critera DOS value to search band edges, cm^-3
# Interpolation type
dos.interp_type = 'linear'
#dos.interp_type = 'cubic'
# integration parameters
# Integration range w.r.t. kBT
dos.Einteg0 = -3.0
dos.Einteg1 = 3.0
dos.nrange = 6.0
# Relative accuracy of the quad functin
dos.eps_quad = 1.0e-5
dos.nmaxdiv_quad = 150
# Bisection method parameters
# Initial search range
dos.dE_bisec = 0.4
# EPS
dos.eps_bisec = 1.0e-7
# max iteration number
dos.nmaxiter_bisec = 200
# Output cycle
dos.iprintiterval_bisec = 1
E_raw = None
dos_raw = None
nDOS = None
Emin = None
Emax = None
Estep = None
EV = None
EC = None
#=============================
# Treat argments
#=============================
[ドキュメント]
def usage():
"""
コマンドライン引数の使用方法を表示する。
詳細説明:
各実行モード（me, T, EF）における引数の指定順序や意味、使用例を標準出力に表示する。
:returns: None
"""
global mode, dos
global Vcell, T0
global file
global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
global dE_me_lsq
print("")
print("Usage: Variables in () are optional")
print(" python {} mode (file args)".format(sys.argv[0], ))
print(" data_for_bandedges: Data to determine band edges [EIGENVAL|DOSCAR]")
print(" (i) python {} me file dE data_for_bandedges".format(sys.argv[0]))
print(" Plot DOS and fit for me*(DOS)")
print(" dE: Energy range from EC/EV for fitting to m*(DOS)")
print(" ex: python {} me {} {} {}".format(sys.argv[0], file, dE_me_lsq, dos.data_for_bandedges))
print(" (ii) python {} T file Tmin Tmax nT data_for_bandedges)".format(sys.argv[0]))
print(" Calculate T dependences of ne, nh, ND+, NA-")
print(" Tmin, Tmax, nT: T range and mesh")
print(" ex: python {} T {} {} {} {} {}".format(sys.argv[0], file, Tmin, Tmax, nT, dos.data_for_bandedges))
print(" (iii) python {} EF file nEF T data_for_bandedges)".format(sys.argv[0]))
print(" Calculate EF dependences of ne, nh, ND+, NA-")
print(" nEF: # of EF mesh")
print(" T: Temperature")
print(" ex: python {} EF {} {} {} {}".format(sys.argv[0], file, nEF, T0, dos.data_for_bandedges))
[ドキュメント]
def updatevars():
"""
コマンドライン引数を解析し、グローバル変数を更新する。
詳細説明:
実行モード、入力ファイル名、および各モードに固有のパラメータを取得し、
出力用CSVファイル名の設定などを行う。
:returns: None
"""
global mode, dos
global Vcell, T0
global file, outcsvfile
global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
global Tmin, Tmax, nT, Tstep
global dE_me_lsq
argv = sys.argv
if len(argv) == 1:
usage()
exit()
mode = getarg( 1, mode)
file = getarg( 2, file)
if mode == 'me':
dE_me_lsq = getfloatarg(3, dE_me_lsq)
dos.data_for_bandedges = getarg(4, dos.data_for_bandedges)
elif mode == 'EF':
nEF = getintarg(3, nEF)
T0 = getfloatarg(4, T0)
dos.data_for_bandedges = getarg(5, dos.data_for_bandedges)
elif mode == 'T':
Tmin = getintarg(3, Tmin)
Tmax = getintarg(4, Tmax)
nT = getintarg(5, nT)
Tstep = (Tmax - Tmin) / (nT - 1)
dos.data_for_bandedges = getarg(6, dos.data_for_bandedges)
header, ext = os.path.splitext(file)
filebody = os.path.basename(header)
outcsvfile = filebody + f'-EFT-{mode}.csv'
[ドキュメント]
def exec_T():
"""
キャリア密度やフェルミ準位の温度(T)依存性を計算してプロットする。
詳細説明:
指定された温度範囲でDOSCARデータに基づき数値積分を行い、フェルミ準位、
電子・正孔密度、不純物密度、活性化エネルギーを算出する。
計算結果はCSVファイルに保存され、Matplotlibを用いてグラフ描画される。
:returns: None
"""
global mode, dos, T0
print("")
print("T dependence")
print(" Temperature range: {} - {}, {} K step".format(Tmin, Tmax, Tstep))
xT = [Tmin + iT * Tstep for iT in range(nT)]
xTinv, yEF, yNe, yNh, yNAm, yNDp, yRH, yNs, yNsabs = dos.cal_T_lists(xT, print_level = 0)
print("")
yEaNe = []
yEaNeTh = []
yEaNh = []
yEaNhTh = []
print ("%8s %10s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s"
% ("T(K)", "EF(eV)", "Ea(Ne,meV)", "Ea(Ne,T1/2)", "Ea(Nh,meV)", "Ea(Nh,T1/2)", "Ne", "Nh", "ND+", "NA-", "RH", "Ns"))
for i in range(nT):
if i == 0:
i1 = i
i2 = i+1
elif i == nT - 1:
i1 = i-1
i2 = i
else:
i1 = i-1
i2 = i+1
Th1 = sqrt(xT[i1])
Th2 = sqrt(xT[i2])
dT = 1.0 / (xTinv[i2] - xTinv[i1])
eps = 1.0e-300
EaNe = (log(yNe[i2]+eps) - log(yNe[i1])+eps) * dT # meV
EaNeTh = (log(yNe[i2]/Th2+eps) - log(yNe[i1]/Th1)+eps) * dT # meV
EaNh = (log(yNh[i2]+eps) - log(yNh[i1])+eps) * dT # meV
EaNhTh = (log(yNh[i2]/Th2+eps) - log(yNh[i1]/Th1)+eps) * dT # meV
yEaNe.append (-1000.0 * 1000.0 * kB / e * EaNe)
yEaNeTh.append(-1000.0 * 1000.0 * kB / e * EaNeTh)
yEaNh.append (-1000.0 * 1000.0 * kB / e * EaNh)
yEaNhTh.append(-1000.0 * 1000.0 * kB / e * EaNhTh)
print("%8.3f %10.4f %12.4g %12.4g %12.4g %12.4g " % (xT[i], yEF[i], yEaNe[i], yEaNeTh[i], yEaNh[i], yEaNhTh[i])
+ "%12g %12g %12g %12g %12g %12g" % (yNe[i], yNh[i], yNDp[i], yNAm[i], yRH[i], yNs[i]))
print("")
print("Save fitting data to [{}]".format(outcsvfile))
save_csv(outcsvfile,
["T(K)", "EF(eV)", "Ne", "Nh", "ND+", "NA-", "RH", "Ns", "Ea(Ne,meV)", "Ea(Ne,T1/2)", "Ea(Nh,meV)", "Ea(Nh,T1/2)"],
[xT, yEF, yNe, yNh, yNDp, yNAm, yRH, yNs, yEaNe, yEaNeTh, yEaNh, yEaNhTh])
#=============================
# Plot graphs
#=============================
fig = plt.figure(figsize = (12, 8))
ax1 = fig.add_subplot(2, 3, 1)
ax2 = fig.add_subplot(2, 3, 2)
ax3 = fig.add_subplot(2, 3, 4)
ax4 = fig.add_subplot(2, 3, 3)
ax5 = fig.add_subplot(2, 3, 6)
ax1.plot(xT, yEF, label = '$E_F$', color = 'black')
ax1.set_xlabel("T (K)")
ax1.set_ylabel("$E_F$ (eV)")
ax1.legend()
ax2.plot(dos.Elist, dos.DOSlist, label = 'DOS', color = 'black')
ax2.plot([dos.EF0, dos.EF0], [min(dos.DOSlist), max(dos.DOSlist)], color = 'red', linestyle = 'dashed')
ax2.set_xlabel("E (eV)")
ax2.set_ylabel("DOS (states/cm$^3$)")
ax2.legend()
ax3.plot(xT, yRH, label = 'RH', color = 'black', marker = 'o')
ax3.set_xlabel("T (K)")
ax3.set_ylabel("R$_H$ (C$^{-1}$cm$^{-3}$)")
ax3.legend()
ax4.plot(xTinv, yNe, label = '$N_e$', color = 'black', marker = 'o', linewidth = 0.5, markersize = 5.0)
ax4.plot(xTinv, yNh, label = '$N_h$', color = 'red', marker = 'o', linewidth = 0.5, markersize = 5.0)
ax4.plot(xTinv, yNsabs, label = '|$N_s$|', color = 'purple', linewidth = 0.5, linestyle = 'dashed', marker = '^', markersize = 2.0)
ax4.plot(xTinv, yNDp, label = '$N_D^+$', color = 'blue', marker = 'x', linewidth = 0.5, markersize = 5.0)
ax4.plot(xTinv, yNAm, label = '$N_A^-$', color = 'green', marker = 'x', linewidth = 0.5, markersize = 5.0)
ax4.set_yscale('log')
minNs = [min(yNe), min(yNh), min(yNsabs)]
maxNs = [max(yNe), max(yNh), max(yNsabs)]
ylim = [ max([min(minNs)*0.5, 1.0e8]), min([max(maxNs)*2.0, 1.0e23]) ]
ax4.set_ylim(ylim)
# ax4.plot([0, 0], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax4.set_xlabel("1000/T (K$^{-1}$)")
ax4.set_ylabel("N (cm$^{-3}$)")
ax4.legend()
ax5.plot(xTinv, yEaNe, label = '$E_a(N_e)$', linewidth = 0.5, color = 'black', marker = 'o', markersize = 5.0)
ax5.plot(xTinv, yEaNeTh, label = '$E_a(N_e*T^{-0.5}$)', linewidth = 0.5, color = 'black', marker = 'x', markersize = 5.0)
ax5.plot(xTinv, yEaNh, label = '$E_a(N_h)$', linewidth = 0.5, color = 'red', marker = 'o', markersize = 5.0)
ax5.plot(xTinv, yEaNhTh, label = '$E_a(N_h*T^{-0.5}$)', linewidth = 0.5, color = 'red', marker = 'x', markersize = 5.0)
ax5.set_xlabel("1000/T (K$^{-1}$)")
ax5.set_ylabel("Ea (meV)")
ax5.legend()
# Rearange the graph axes so that they are not overlapped
plt.tight_layout()
plt.pause(0.1)
print("")
usage()
print("")
print("Press ENTER to exit>>", end = '')
input()
[ドキュメント]
def exec_EF():
"""
キャリア密度等のフェルミ準位(EF)依存性を計算してプロットする。
詳細説明:
特定の温度においてフェルミ準位を変化させ、電子・正孔密度や有効状態密度(DOS)を計算する。
自由電子近似に基づく近似式との比較を行い、結果をCSV保存およびグラフ描画する。
:returns: None
"""
global mode, dos, T0
global EFmin, EFmax, EFstep
print("")
print("EF dependence at {} K".format(T0))
print(" EF range: {} - {}, {} eV step".format(EFmin, EFmax, EFstep))
xEF = [EFmin + iEF * EFstep for iEF in range(nEF)]
yNh, yNe, yNAm, yNDp, yRH, yNs, yNsabs = dos.cal_EF_lists(T0, xEF, print_level = 0)
yT0Ne = []
yT0Nh = []
print("")
print ("%5s\t%8s\t%12s\t%12s\t%12s\t%12s\t%12s\t%12s\t%12s" % ("EF(eV)", "Ne", "Nh", "ND+", "NA-", "RH", "T0(Ne)", "T0(Nh)", "Ns"))
for i in range(nEF):
EF = xEF[i]
if i == 0:
i1 = i
i2 = i+1
elif i == nEF - 1:
i1 = i-1
i2 = i
else:
i1 = i-1
i2 = i+1
dEF = xEF[i2] - xEF[i1]
dlnNedEF = (log(yNe[i2]+1.0e-300) - log(yNe[i1]+1.0e-300)) / dEF
dlnNhdEF = (log(yNh[i2]+1.0e-300) - log(yNh[i1]+1.0e-300)) / dEF
if dlnNedEF == 0.0:
dlnNedEF = 1.0-4
T0Ne = 1.0 / kB / dlnNedEF * e
if dlnNhdEF == 0.0:
dlnNhdEF = 1.0-4
T0Nh = -1.0 / kB / dlnNhdEF * e
yT0Ne.append(T0Ne)
yT0Nh.append(T0Nh)
print ("%10.4f\t%10.4g\t%10.4g\t%10.4g\t%10.4g\t%10.4g\t%10.4g\t%10.4g\t%10.4g"
% (EF, yNe[i], yNh[i], yNDp[i], yNAm[i], yRH[i], T0Ne, T0Nh, yNs[i]))
Emid = (dos.EV + dos.EC) / 2.0
imid = int((Emid - EFmin) / EFstep + 1.0e-5)
Eimid = EFmin + imid * EFstep
Ne = yNe[imid]
Nh = yNh[imid]
T0Ne = yT0Ne[imid]
T0Nh = yT0Ne[imid]
NC = Ne * exp((dos.EC - Eimid) * e / kB / T0Ne)
NV = Nh * exp((Eimid - dos.EV) * e / kB / T0Nh)
meeff = NC2meff_FEA(NC, T0)
mheff = NC2meff_FEA(NV, T0)
DC0 = meff2DC0_FEA(meeff, T0)
DV0 = meff2DC0_FEA(mheff, T0)
DOS_FEA = []
print("")
print ("{:10}\t{:10}".format("E(eV)", "DOS,FEA(1/cm^3/eV)"))
for i in range(nEF):
EF = xEF[i]
if EF < dos.EV:
DOS_FEA.append(DV0 * sqrt(dos.EV - EF))
elif EF > dos.EC:
DOS_FEA.append(DC0 * sqrt(EF - dos.EC))
else:
DOS_FEA.append(0.0)
print ("{:10.4g}\t{:10.4g}".format(xEF[i], DOS_FEA[i]))
print("")
print("Effective density of states at the mid gap {} eV (i={}, E={:8.4g}):".format(Emid, imid, Eimid))
print(" NC = {:12.6g} cm-3 (T0 = {:8.4g} K)".format(NC, T0Ne))
print(" NV = {:12.6g} cm-3 (T0 = {:8.4g} K)".format(NV, T0Nh))
print(" DC0 = {:12.6g} cm^-3/eV^1.5 (T0 = {:8.4g} K)".format(DC0, T0Ne))
print(" DV0 = {:12.6g} cm^-3/eV^1.5 (T0 = {:8.4g} K)".format(DV0, T0Nh))
print(" me* = {:8.4g} me".format(meeff))
print(" mh* = {:8.4g} me".format(mheff))
xEFapprox = [EFmin + i * EFstep for i in range(nEF)]
yNeapprox = [NC * exp(-(dos.EC - xEFapprox[i]) * e / kB / T0) for i in range(nEF)]
yNhapprox = [NV * exp(-(xEFapprox[i] - dos.EV) * e / kB / T0) for i in range(nEF)]
print("")
print("Save fitting data to [{}]".format(outcsvfile))
save_csv(outcsvfile,
["EF(eV)", "Ne", "Nh", "ND+", "NA-", "RH", "T0(Ne)", "T0(Nh)", "Ns"],
[xEF, yNe, yNh, yNDp, yNAm, yRH, yT0Ne, yT0Nh, yNs])
#=============================
# Plot graphs
#=============================
fig = plt.figure(figsize = (8, 8))
ax2 = fig.add_subplot(2, 2, 1)
ax3 = fig.add_subplot(2, 2, 2)
ax4 = fig.add_subplot(2, 2, 3)
ax1 = fig.add_subplot(2, 2, 4)
ax2.plot(dos.Elist, dos.DOSlist, label = 'DOS', color = 'black')
ax2.plot(xEF, DOS_FEA, label = 'DOS$_{,FEA}$', color = 'green', linewidth = 0.5, linestyle = 'dashed')
ax2.plot([dos.EF0, dos.EF0], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax2.plot([dos.EV, dos.EV], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax2.plot([dos.EC, dos.EC], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax2.set_xlabel("E (eV)")
ax2.set_ylabel("DOS (states/cm$^3$)")
ax2.legend()
ax3.plot(xEF, yNe, label = '$N_e$', color = 'black', linewidth = 0.5) #, marker = 'o')
ax3.plot(xEF, yNh, label = '$N_h$', color = 'red', linewidth = 0.5) #, marker = 'x')
ax3.plot(xEF, yNsabs, label = '|$N_s$|', color = 'purple', linewidth = 0.5, linestyle = 'dashed') #, marker = '^')
ax3.plot(xEF, yNDp, label = '$N_D^+$', color = 'blue') #, marker = 'x')
ax3.plot(xEF, yNAm, label = '$N_A^-$', color = 'green') #, marker = 'x')
ax3.plot(xEFapprox, yNeapprox, label = 'Ne(Boltz)', color = 'gray', linewidth = 0.5) #, marker = 'x')
ax3.plot(xEFapprox, yNhapprox, label = 'Nh(Boltz)', color = 'purple', linewidth = 0.5) #, marker = 'x')
ax3.plot([dos.EV, dos.EV], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax3.plot([dos.EC, dos.EC], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax3.set_yscale('log')
ax3.set_xlabel("EF (eV)")
ax3.set_ylabel("N (cm$^{-3}$)")
ax3.legend()
ax4.plot(xEF, yRH, label = 'R$_H$', color = 'black', linewidth = 0.5, marker = 'o')
ax4.set_xlabel("EF (eV)")
ax4.set_ylabel("R$_H$ (C$^{-1}$cm$^{-3}$)")
ax4.legend()
ax1.plot(xEF, yT0Ne, label = 'T$_0$(Ne)', color = 'black') #, marker = 'o')
ax1.plot(xEF, yT0Nh, label = 'T$_0$(Nh)', color = 'red') #, marker = 'x')
ax1.plot([dos.EV, dos.EV], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax1.plot([dos.EC, dos.EC], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax1.set_ylim([0.0, 2.0*T0])
ax1.set_xlabel("EF (eV)")
ax1.set_ylabel("T$_0$ (K)")
ax1.legend()
# Rearange the graph axes so that they are not overlapped
plt.tight_layout()
plt.pause(0.1)
print("")
usage()
print("")
print("Press ENTER to exit>>", end = '')
input()
[ドキュメント]
def exec_me():
"""
状態密度(DOS)をプロットし、バンド端付近の有効質量(m*)を見積もる。
詳細説明:
DOSの2乗とエネルギーの関係から傾きを計算し、自由電子近似に基づく
有効質量を算出する。各エネルギーごとのDOSやm*の依存性をプロットする。
:returns: None
"""
global dos
global file
print("")
print("Plot DOS and calculate m*")
print(" dE_me_lsq: {:8.4g} eV".format(dE_me_lsq))
print(" EV: {:8.4g} eV dEFmin: {:8.4g} eV".format(dos.EV, dEFmin))
print(" EC: {:8.4g} eV dEFmax: {:8.4g} eV".format(dos.EC, dEFmax))
plotE0 = min([dos.EV + dEFmin, dos.EV - dE_me_lsq])
plotE1 = max([dos.EC + dEFmax, dos.EC + dE_me_lsq])
print("Plot E range: {:8.4g} - {:8.4g} eV".format(plotE0, plotE1))
ne_from_dos = integrate_Simpson_list(dos.Elist, dos.DOSlist)
E_me = []
dos2 = []
Emid = (dos.EV + dos.EC) / 2.0
c = 0
for i in range(dos.nE):
E = dos.Elist[i]
# print("ne=", dos.Emin, E, dos.DOS(E), ne)
if E < plotE0 or plotE1 < E:
continue
E_me.append(E)
dos2.append(dos.DOSlist[i]**2)
me_me = []
nme = len(E_me)
for i in range(nme):
if i == 0:
i0 = 0
i1 = 1
elif i == nme - 1:
i0 = i - 1
i1 = i
else:
i0 = i - 1
i1 = i + 1
DC0 = (dos2[i1] - dos2[i0]) / (E_me[i1] - E_me[i0])
if DC0 < 0.0:
sign = -1.0
else:
sign = 1.0
DC0 = sqrt(abs(DC0))
# DC0 = sqrt(2.0) / pi / pi * pow(me * meff, 1.5) / hbar / hbar / hbar * 1.0e-6 * pow(e, 1.5)
meff = pow(DC0 / sqrt(2.0) * pi * pi * hbar * hbar * hbar / 1.0e-6 / pow(e, 1.5), 2.0/3.0) / me
meff *= sign
# print("m=", DC0, meff)
me_me.append(meff)
#=============================
# Plot graphs
#=============================
fig = plt.figure(figsize = (8, 8))
ax1 = fig.add_subplot(2, 2, 1)
ax2 = fig.add_subplot(2, 2, 2)
ax3 = fig.add_subplot(2, 2, 3)
ax4 = fig.add_subplot(2, 2, 4)
ax1.plot(dos.Elist, dos.DOSlist, label = 'DOS', color = 'black', linewidth = 0.5)
ax1.plot([dos.EF0, dos.EF0], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax1.plot([dos.EV, dos.EV], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax1.plot([dos.EC, dos.EC], [min(dos.DOSlist), max(dos.DOSlist)], label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax1.set_xlabel("E (eV)")
ax1.set_ylabel("DOS (states/cm$^3$/eV)")
ax1.legend()
# print("len=", len(dos.Elist), len(ne_from_dos))
# print("ne=", ne_from_dos)
ax2.plot(dos.Elist, ne_from_dos, label = 'N', color = 'black', linewidth = 0.5)
ylim = [min(ne_from_dos), max(ne_from_dos)]
ax2.plot([dos.EF0, dos.EF0], ylim, label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax2.plot([dos.EV, dos.EV], ylim, label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax2.plot([dos.EC, dos.EC], ylim, label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax2.set_xlabel("E (eV)")
ax2.set_ylabel("N (states/cm$^3$)")
ax2.legend()
ax3.plot(E_me, dos2, label = 'DOS$^2$', color = 'black', linewidth = 0.5)
ylim = [min(dos2), max(dos2)]
ax3.plot([dos.EF0, dos.EF0], ylim, label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax3.plot([dos.EV, dos.EV], ylim, label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax3.plot([dos.EC, dos.EC], ylim, label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
# ax3.plot([dos.EV - dE_me_lsq, dos.EV - dE_me_lsq], ylim, label = '$E_{fit}$', color = 'green', linestyle = 'dashed', linewidth = 0.5)
# ax3.plot([dos.EC + dE_me_lsq, dos.EC + dE_me_lsq], ylim, label = '$E_{fit}$', color = 'green', linestyle = 'dashed', linewidth = 0.5)
ax3.set_xlabel("E (eV)")
ax3.set_ylabel("DOS$^2$ (states$^2$/cm$^6$)")
ax3.legend()
ax4.plot(E_me, me_me, label = '$m^*$', color = 'black', linewidth = 0.5)
ylim = [min(me_me), max(me_me)]
ax4.plot([dos.EF0, dos.EF0], ylim, label = '$E_{F0}$', color = 'red', linestyle = 'dashed', linewidth = 0.5)
ax4.plot([dos.EV, dos.EV], ylim, label = '$E_V$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
ax4.plot([dos.EC, dos.EC], ylim, label = '$E_C$', color = 'blue', linestyle = 'dashed', linewidth = 0.5)
# ax4.plot([dos.EV - dE_me_lsq, dos.EV - dE_me_lsq], ylim, label = '$E_{fit}$', color = 'green', linestyle = 'dashed', linewidth = 0.5)
# ax4.plot([dos.EC + dE_me_lsq, dos.EC + dE_me_lsq], ylim, label = '$E_{fit}$', color = 'green', linestyle = 'dashed', linewidth = 0.5)
ax4.set_xlabel("E (eV)")
ax4.set_ylabel("$m^*$")
ax4.legend()
# Rearange the graph axes so that they are not overlapped
plt.tight_layout()
plt.pause(0.1)
print("")
usage()
print("")
print("Press ENTER to exit>>", end = '')
input()
[ドキュメント]
def main():
"""
プログラムのメインルーチンを実行する。
詳細説明:
変数の初期化と更新を行い、VASPのDOSCARファイルを読み込んだ後、
指定された実行モード（me, T, EF）に対応する処理を呼び出す。
:returns: None
"""
global mode, dos, T0
global Vcell, file, E_raw, dos_raw
global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
updatevars()
print("")
print("===============================================================================")
print(" Calculate T dependences of semiconductor statistics using VASP results")
print("===============================================================================")
print("mode: ", mode)
dos.read_vaspfiles(file)
EFmin = dos.EV + dEFmin
EFmax = dos.EC + dEFmax
EFstep = (EFmax - EFmin) / (nEF - 1)
print("")
print("Integration configuration")
print(" Integration E range: {} - {} eV, or EF+-{}*kBT eV".format(dos.Einteg0, dos.Einteg1, dos.nrange))
print(" quad() parameters")
print(" releps=", dos.eps_quad)
print(" nmaxdiv=", dos.nmaxdiv_quad)
print("")
print("EF at 0 K = ", dos.EF0, " eV")
Ne0K = dos.Ne(0.0, dos.EF0, -5.0, dos.EF0, dos.Estep)
# Ne0K = Ne(0.0, EF0, Einteg0, EF0)
print(" Ne at 0 K = ", Ne0K, " /u.c.")
print("")
if mode == 'me':
exec_me()
elif mode == 'T':
exec_T()
elif mode == 'EF':
exec_EF()
else:
terminate("Error: Invalid mode [{}]".format(mode), usage)
if __name__ == '__main__':
main()