"""
バンド構造データから有効質量を計算しプロットするスクリプトです。
このスクリプトは、k空間でのエネルギーバンドデータ（CSV形式）を読み込み、
その二次微分を用いて有効質量を計算します。計算された有効質量は、
自由電子質量meに対する相対値としてプロットされます。
特に、有効質量が極端に大きくなる（逆有効質量が0に近づく）点を
識別し、プロット上で区切りとして表現します。
関連リンク:
:doc:`EffectiveMass_usage`
"""
import sys
import numpy as np
from numpy import sqrt, exp, sin, cos, tan, pi
import numpy.linalg as LA
import csv
from matplotlib import pyplot as plt
#===================================
# physical constants
#===================================
pi = 3.14159265358979323846
pi2 = 2.0 * pi
h = 6.6260755e-34 # Js";
hbar = 1.05459e-34 # "Js";
c = 2.99792458e8 # m/s";
e = 1.60218e-19 # C";
e0 = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
kB = 1.380658e-23 # JK<sup>-1</sup>";
me = 9.1093897e-31 # kg";
R = 8.314462618 # J/K/mol
a0 = 5.29177e-11 # m";
#===================================
# parameters
#===================================
a = 4.0e-10 #m
infile = 'band.csv'
cutline = 1
#===================================
# figure configuration
#===================================
fontsize = 12
legend_fontsize = 8
#=============================
# other functions
#=============================
[ドキュメント]
def savecsv(outfile, header, datalist):
"""
指定されたヘッダーとデータリストをCSVファイルに保存します。
`datalist`は列指向のデータ構造として扱われ、各内部リストが1つの列を表します。
ファイルへの書き込み時には、これらの列が転置されて行として書き込まれます。
ファイル書き込みに失敗した場合はエラーメッセージを出力します。
:param outfile: str: 出力するCSVファイルのパス。
:param header: list[str]: CSVファイルのヘッダー行として使用される文字列のリスト。
:param datalist: list[list[float]]: 保存するデータ。
datalist[j][i]がi行j列の要素となるような（列優先の）構造を想定しています。
例: `[[col1_row1, col1_row2], [col2_row1, col2_row2]]`
:returns: None
"""
try:
print("Write to [{}]".format(outfile))
f = open(outfile, 'w')
except:
# except IOError:
print("Error: Can not write to [{}]".format(outfile))
else:
fout = csv.writer(f, lineterminator='\n')
fout.writerow(header)
# fout.writerows(data)
for i in range(0, len(datalist[0])):
a = []
for j in range(len(datalist)):
a.append(datalist[j][i])
fout.writerow(a)
f.close()
[ドキュメント]
def read_csv(fname):
"""
2列のCSVファイルを読み込み、x軸/y軸のラベルとデータを返します。
最初の行をヘッダーとしてx軸とy軸のラベルを抽出し、それ以降の行を
浮動小数点数データとして読み込みます。数値変換に失敗した行は警告を
出力し、そのデータはスキップされます。
:param fname: str: 読み込むCSVファイルのパス。
:returns: tuple[str, str, list[float], list[float]]:
x軸のラベル、y軸のラベル、xデータ（浮動小数点数のリスト）、yデータ（浮動小数点数のリスト）。
"""
x = []
y = []
with open(fname) as f:
fin = csv.reader(f)
aa = next(fin)
xlabel = aa[0]
ylabel = aa[1]
for row in fin:
try:
x.append(float(row[0]))
y.append(float(row[1]))
except:
print("Warning: Invalid float data [{}] or [{}]".format(row[0], row[1]))
return xlabel, ylabel, x, y
[ドキュメント]
def main():
"""
メイン処理を実行し、有効質量の計算とプロットを行います。
`infile`で指定されたCSVファイルからk-Eバンド構造データを読み込み、
中心差分法を用いてエネルギーEのkに関する2次微分d^2E/dk^2を計算します。
この2次微分から有効質量の逆数(1/m*)を導出し、自由電子質量meに対する
相対有効質量(m*/me)を求めます。
有効質量が極端に大きくなる（逆有効質量がほぼ0になる）点や、有効質量の
符号が反転する点では、プロットを一時的に切断し、視認性を高めます。
`nskip`の異なる2つのケース（`nskip=1`と`nskip=4`）で計算とプロットを行い、
計算結果を比較できるようにします。最後に、結果をMatplotlibで表示します。
:returns: None
"""
klabel, Elabel, k, E = read_csv(infile)
print("k=", k)
print("E=", E)
nk = len(k)
dk = k[1] - k[0]
km = hbar * hbar * (pi2 / a)**2.0
nskip = 1
xk = []
ymc = []
signprev = None
print("")
print("Differentiate for nskip = ", nskip)
print("i E[i-1] E[i] E[i+1] 1.0/m*")
for i in range(nskip, nk - nskip, nskip):
# 2次微分 d^2E/dk^2 の計算 (中心差分法)
# 物理量に変換するため電荷eを掛ける
d2Edk2c = (E[i+nskip] + E[i-nskip] - 2 * E[i]) * e / pow(nskip * dk, 2.0)
minv = d2Edk2c / km # 逆有効質量
print(i, E[i-1], E[i], E[i+1], minv)
# 逆有効質量が極端に小さい場合 (有効質量が極端に大きい場合) は、プロットを切断する
if abs(minv) <= 1.0e20: # << 1.0/me ~ 1e30
if cutline:
xk.append(k[i])
ymc.append(None) # Noneを挿入してプロットを分断
signprev = -signprev # 符号反転を強制
continue
else:
m = km / d2Edk2c # 有効質量を計算
# 有効質量の符号が反転した場合もプロットを切断する
if signprev is None:
signprev = m
elif signprev * m < 0.0:
if cutline:
xk.append(k[i])
ymc.append(None) # Noneを挿入してプロットを分断
signprev = m
xk.append(k[i])
ymc.append(m / me) # 自由電子質量で規格化して追加
nskip = 4
xk2 = []
ymc2 = []
signprev = None
print("")
print("Differentiate for nskip = ", nskip)
print("i E[i-1] E[i] E[i+1] 1.0/m*")
for i in range(nskip, nk - nskip, nskip):
d2Edk2c = (E[i+nskip] + E[i-nskip] - 2 * E[i]) * e / pow(nskip * dk, 2.0)
minv = d2Edk2c / km
print(i, E[i-1], E[i], E[i+1], minv)
if abs(minv) <= 1.0e20: # << 1.0/me ~ 1e30
if cutline:
xk2.append(k[i])
ymc2.append(None)
signprev = -signprev
continue
else:
m = km / d2Edk2c
if signprev is None:
signprev = m
elif signprev * m < 0.0:
if cutline:
xk2.append(k[i])
ymc2.append(None)
signprev = m
xk2.append(k[i])
ymc2.append(m / me)
# 結果のプロット
plt.plot(xk, ymc, linewidth = 0.5, marker = 'o', markersize = 1.0, label = 'nskip = 1')
plt.plot(xk2, ymc2, linewidth = 0.5, marker = 'o', markersize = 1.0, label = 'nskip = 2') # ラベルが"nskip = 2"だが、コードではnskip=4で計算しているため注意
plt.xlabel(klabel)
plt.ylabel("m$_e$ / m$_e^0$")
plt.xlim([-0.5, 0.5])
# plt.ylim([-0.5, 0.5])
plt.tight_layout()
plt.legend(fontsize=legend_fontsize) # 凡例を表示
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
if __name__ == "__main__":
main()