問題
バンド構造 band.csv から、数値微分により d2E(k)/dk2 を求め、有効質量 me*と k の関係をグラフに描け。
格子定数は a = 4.0 Å とする。
異なる精度の数値微分をし、有効質量の精度の比較をするとbetter。

有効質量TBバンドE (k ) = ε 1 − 2 h12 cos(ka ) ~ ε 1 − 2 h12 + h12 a 2 k 2 + O((ka ) )42自由電子P
2 2= E0 +E (k ) = E0 +k2m2m

11 ∂ 2 E n (k )= 2*m∂k 2

Energy / |h12|22m* =
2 h12 a 2

大きな混成(|h12|)により質量m*は小さくなる

2m = 2 / Wa

10.50

バンド幅 W = 4|h12|*e1.5-0.5

2-0.3-0.1

2π/a0.10.30.5有効質量k は逆格子の内部座標: 一般に [-½ ½] の範囲で表示される
単位変換 kreal = (2π/a) k−1 ∂ En (k ) 2  2π   ∂ En (k )  = 
m =   2 2
 a   ∂k ∂k real *2222−1

一般に、有効質量は 電子の静止質量 me0 との比であらわす−1−1 ∂ En (k ) 002  2π   ∂ En (k ) m / me =  /m/m= ee2 2ak∂  ∂k real *02222

数値計算: 微分df ( x)をコンピュータでどのように計算するか
dx

微分 d を差分 Δ で置き換える
df ( x) ∆f ( x) f ( x + h) − f ( x) f ( x + h) − f ( x)
~
=
=
∆xhdx
( x + h) − xh を 小さくすれば精度が上がる  桁落ち誤差
32bit浮動小数点 (~7桁) : 扱う最小数値の 5桁下が限界
64bit浮動小数点 (~16桁): 扱う最小数値の 14桁下が限界
df ( x)1 d 2 f ( x) 2
3f ( x + h) = f ( x ) +h+h+Oh()
2dx
2 dx
f ( x + h) − f ( x) df ( x) 1 d 2 f ( x)2=+h+Oh() 差分誤差
2hdx
2 dx

数値微分: 平均を取って精度を上げる
df ( x) f ( x + h) − f ( x)~hdx
1 d 3 f ( x) 2f ( x + h) − f ( x) df ( x) 1 d 2 f ( x)4()=+h+h+Oh
誤差：
232 dx
3! dx
hdx
f ( x + h) − f ( x − h)df ( x)  f ( x + h) − f ( x) f ( x) − f ( x − h) ~+/ 2=
2hhhdx
df ( x)1 d 2 f ( x) 2 1 d 3 f ( x) 34f ( x + h) = f ( x ) +h+hh
Oh()++
23dx
2 dx
3! dx
df ( x)1 d 2 f ( x) 2 1 d 3 f ( x) 3
4f ( x − h) = f ( x ) −h+hhOh()−+23dx2 dx
3! dxf ( x + h) − f ( x − h) df ( x) 1 d 3 f ( x) 34h
Oh()
=++
誤差：
3dx
2h
3! dx

二階微分
一階微分を前進差分で計算してから二階微分を計算すると・・・

dx
dx
()
(t )
tt+
∆−d 2 x (t ) dtdt=dt 2∆tx(t + ∆t ) − x(t ) x(t ) − x(t − ∆t )−x(t + 2∆t ) − 2 x(t + ∆t ) + x(t )
∆t
∆t~=
∆t
∆t 2
偶数階微分では、結果が t + Δt、t – Δt について対称になる式を取れる

x(t + ∆t ) − x(t ) x(t ) − x(t − ∆t )−d 2 x (t )∆t∆t~2dt∆tx(t + ∆t ) − 2 x(t ) + x(t − ∆t )=∆t 2

２つの式では、横軸がひとつずれるので注意！

プログラム (抜粋)

EffectiveMass.py

描いてはいけない線

青: nskip =1
橙: nskip =4

データに None (未定義値) を挿入することで描いてはいけない線を消した

青: nskip =1
橙: nskip =4

EffectiveMass.py

プログラム (抜粋)#===================================# parameters#===================================#a = 4.0 # A
a = 4.0e-10 #m#共通の定数は先に計算
km = hbar * hbar * (pi2 / a)**2.0#微分の精度を比較するため、h = nskip*dk にする
nskip = 1
xk = []
ymc = []
infile = 'band.csv'#符号の変化を検出するため、符号変数を用意
signprev = None
#有効質量の符号が変わる点をつなぐかどうかのフラグ
for i in range(nskip, nk - nskip, nskip):cutline = 1#2階微分を計算
d2Edk2c = (E[i+nskip] + E[i-nskip] - 2 * E[i]) * e / pow(nskip *def read_csv(fname):x = []
dk, 2.0)y = []#2回微分はゼロになることがあるので、まずは1/m*を計算
minv = d2Edk2c / kmwith open(fname) as f:fin = csv.reader(f)print(i, E[i-1], E[i], E[i+1], minv)
xlabel, ylabel, = next(fin)#1/m*が1/meより非常に小さければ、m*は計算しない
for row in fin:
if abs(minv) <= 1.0e20: # << 1.0/me ~ 1e30
try:
#符号が反転する場所でグラフの線を切断するときは
x.append(float(row[0]))#Noneデータを追加する。
y.append(float(row[1]))if cutline:
except:xk.append(k[i])print("Warning: Invalid float data [{}] or [{}]".format(row[0],ymc.append(None)row[1]))#反転した符号を記録
signprev = -signprev
return xlabel, ylabel, x, y
continueelse:
def main():m = km / d2Edk2c
klabel, Elabel, k, E = read_csv(infile)
if signprev is None:#入力データで使う変数を計算
nk = len(k)signprev = mdk = k[1] - k[0]elif signprev * m < 0.0:
if cutline:EffectiveMass.py

プログラム (抜粋)

#共通の定数は先に計算
km = hbar * hbar * (pi2 / a)**2.0#微分の精度を比較するため、h = nskip*dk にする
nskip = 1
xk = []
ymc = []
#符号の変化を検出するため、符号変数を用意
signprev = Nonefor i in range(nskip, nk - nskip, nskip):#2階微分を計算
d2Edk2c = (E[i+nskip] + E[i-nskip] - 2 * E[i]) * e / pow(nskip *dk, 2.0)
#2回微分はゼロになることがあるので、まずは1/m*を計算
minv = d2Edk2c / km
print(i, E[i-1], E[i], E[i+1], minv)#1/m*が1/meより非常に小さければ、m*は計算しない
if abs(minv) <= 1.0e20: # << 1.0/me ~ 1e30
#符号が反転する場所でグラフの線を切断するときは#Noneデータを追加する。
if cutline:
xk.append(k[i])ymc.append(None)#反転した符号を記録
signprev = -signprev
continueelse:
m = km / d2Edk2c#符号が反転する場所でグラフの線を切断するときは#Noneデータを追加する。
if signprev is None:#signprevが 初期値 None である場合は 符号の最初の値を代入
signprev = melif signprev * m < 0.0:
if cutline:xk.append(k[i])ymc.append(None)#反転した符号を記録
signprev = m
xk.append(k[i])ymc.append(m / me)plt.plot(xk, ymc, linewidth = 0.5, marker = 'o', markersize = 1.0,
label = 'nskip = 1')
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.pause(0.1)
print("Press ENTER to exit>>", end = '')input()
if __name__ == "__main__":main()