import sys
from pprint import pprint
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
"""
概要: 自由電子バンド構造計算を行うスクリプト。
詳細説明:
このスクリプトは、自由電子近似に基づき物質の電子バンド構造を計算し、プロットします。
結晶の格子定数、有効質量、逆格子ベクトル、計算するk点パスなどを設定し、
高対称点に沿ったバンド分散 E(k) を計算・可視化します。
関連リンク:
:doc:`free_electron_band_usage`
"""
#===================================
# 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";
#========================
# Crystal definition
#========================
# Si
a = 5.4064 # angstrom, lattice parameter
rg = np.zeros([3, 3])
rg[0][0] = pow(2.0 * pi / a, 2) # reciprocal space metric in angstrom^-2
rg[1][1] = rg[0][0]
rg[2][2] = rg[0][0]
# 有効質量
meff = 1.0 # in me
# E(k) = KE * k*k, E(k) in eV, k in angstrom^1
KE = hbar * hbar / 2.0 / (meff * me) * 1.0e20 / e
print("KE = ", KE)
#========================
# Band plot parameters
#========================
# バンド構造をプロットするk点の軌跡: [kx, ky, kz, k点名称]
klist = [
# [-0.5, 0.0, 0.0, "-X"]
# , [0.0, 0.0, 0.0, "$\Gamma$"]
# , [ 0.5, 0.0, 0.0, "X"]
[0.5, 0.0, 1.0, "W"]
, [0.5, 0.5, 0.5, "L"]
, [0.0, 0.0, 0.0, "$\Gamma$"]
, [0.0, 0.0, 1.0, "X"]
, [0.5, 0.0, 1.0, "W"]
, [0.75, 0.0, 0.75, "K"]
]
# プロットするバンド構造E(k)のk点数の概数
nk = 101
# Ehkl(k)を計算するhkl範囲
hrange = [-3, 3]
krange = [-3, 3]
lrange = [-3, 3]
#hrange = [0, 0]
#krange = [0, 0]
#lrange = [0, 0]
# プロットするエネルギー範囲
Erange = [0.0, 10.0] # eV
#===================================
# figure configuration
#===================================
figsize = (6, 8)
#figsize = (6, 4)
fontsize = 16
legend_fontsize = 8
#==============================================
# fundamental functions
#==============================================
[ドキュメント]
def pfloat(str):
"""
概要: 文字列を浮動小数点数に安全に変換します。
詳細説明:
`float()` 変換に失敗した場合でもプログラムを終了させず、Noneを返します。
:param str: str 変換する文字列。
:returns: float or None 変換された浮動小数点数、または変換できなかった場合はNone。
"""
try:
return float(str)
except:
return None
[ドキュメント]
def pint(str):
"""
概要: 文字列を整数に安全に変換します。
詳細説明:
`int()` 変換に失敗した場合でもプログラムを終了させず、Noneを返します。
:param str: str 変換する文字列。
:returns: int or None 変換された整数、または変換できなかった場合はNone。
"""
try:
return int(str)
except:
return None
[ドキュメント]
def getarg(position, defval = None):
"""
概要: コマンドライン引数を安全に取得します。
詳細説明:
指定されたインデックスのコマンドライン引数を取得します。
インデックスが範囲外の場合は、デフォルト値を返します。
:param position: int 取得する引数のインデックス（0はスクリプト名）。
:param defval: Any, optional 取得できなかった場合に返すデフォルト値。デフォルトはNone。
:returns: Any コマンドライン引数の値、またはデフォルト値。
"""
try:
return sys.argv[position]
except:
return defval
[ドキュメント]
def getfloatarg(position, defval = None):
"""
概要: コマンドライン引数を浮動小数点数として安全に取得します。
詳細説明:
`getarg` を使用して引数を取得し、`pfloat` で浮動小数点数に変換します。
:param position: int 取得する引数のインデックス。
:param defval: float or None, optional 取得または変換できなかった場合に返すデフォルト値。デフォルトはNone。
:returns: float or None 変換された浮動小数点数、またはデフォルト値。
"""
return pfloat(getarg(position, defval))
[ドキュメント]
def getintarg(position, defval = None):
"""
概要: コマンドライン引数を整数として安全に取得します。
詳細説明:
`getarg` を使用して引数を取得し、`pint` で整数に変換します。
:param position: int 取得する引数のインデックス。
:param defval: int or None, optional 取得または変換できなかった場合に返すデフォルト値。デフォルトはNone。
:returns: int or None 変換された整数、またはデフォルト値。
"""
return pint(getarg(position, defval))
[ドキュメント]
def usage():
"""
概要: プログラムの利用方法（usage）メッセージを表示します。
詳細説明:
現在ではシンプルなメッセージを表示しますが、将来的にはより詳細な説明を追加できます。
"""
print("")
print("Usage:")
# print(" python {}".format(sys.argv[0]))
[ドキュメント]
def terminate(message = None):
"""
概要: メッセージを表示してプログラムを終了します。
詳細説明:
オプションのメッセージを表示し、`usage` 関数を呼び出した後、プログラムを終了します。
:param message: str or None, optional 表示する終了メッセージ。デフォルトはNone。
"""
print("")
if message is not None:
print("")
print(message)
print("")
usage()
print("")
exit()
[ドキュメント]
def cal_kdistance(rg, k0, k1):
"""
概要: 逆格子空間における2つのk点間の距離を計算します。
詳細説明:
与えられた逆格子のmetrics (rg) を用いて、2つのk点ベクトル間の距離を計算します。
k点ベクトルは内部座標で与えられます。
:param rg: numpy.ndarray 逆格子のmetric tensor (3x3)。
:param k0: list or numpy.ndarray 最初のk点ベクトル [kx, ky, kz]。
:param k1: list or numpy.ndarray 2番目のk点ベクトル [kx, ky, kz]。
:returns: float 2つのk点間の距離。
"""
dkx = k1[0] - k0[0]
dky = k1[1] - k0[1]
dkz = k1[2] - k0[2]
r2 = rg[0][0] * dkx*dkx + rg[1][1] * dky*dky + rg[2][2] * dkz*dkz
r2 += 2.0 * (rg[0][1] * dkx*dky + rg[1][2] * dky*dkz + rg[2][0] * dky*dkx)
return sqrt(r2)
[ドキュメント]
def cal_E(k, Ghkl):
"""
概要: 自由電子モデルにおけるエネルギーE(k)を計算します。
詳細説明:
与えられたk点ベクトルと逆格子ベクトルG (hkl) を用いて、自由電子のエネルギーを計算します。
k点とGhklは内部座標系で与えられ、エネルギーはeV単位で返されます。
:param k: list or numpy.ndarray 内部座標系でのk点ベクトル [kx, ky, kz]。
:param Ghkl: list or numpy.ndarray 逆格子ベクトル [h, k, l] (ミラー指数)。
:returns: float 計算された自由電子のエネルギー (eV)。
"""
global rg
kabs2 = rg[0][0] * (k[0] + Ghkl[0])**2 # in angstrom^-2
kabs2 += rg[1][1] * (k[1] + Ghkl[1])**2
kabs2 += rg[2][2] * (k[2] + Ghkl[2])**2
E = KE * kabs2 # in eV
return E
[ドキュメント]
def get_dklist(klist, nk):
"""
概要: プロットするk点パスにおける、各セグメントのk点距離と累積距離を計算します。
詳細説明:
`klist` に定義された高対称点パスに沿って、各点間の距離 (`dklist`) および
最初のk点からの累積距離 (`ktotal_list`) を計算します。
また、バンド構造プロット用のk点名称リスト (`ktotal_namelist`) も生成します。
:param klist: list k点パスのリスト。各要素は [kx, ky, kz, k点名] の形式。
:param nk: int 計算するk点の概数。この関数では直接使用されないが、バンド計算の粒度を示すため引数として残す。
:returns dklist: list 各k点セグメントの距離のリスト。
:returns ktotal_list: list 最初のk点からの累積距離のリスト。
:returns ktotal_namelist: list 累積距離リストに対応するk点名称のリスト。
:returns ktotal: float 全てのk点パスの総距離。
"""
print("")
# list of k distances from the first k point of each k region
dklist = []
# list of k distance from the first k point of the k list
ktotal = 0.0
ktotal_list = []
ktotal_namelist = []
ktotal_list.append(0.0)
ktotal_namelist.append(klist[0][3])
for i in range(1, len(klist)):
print("k [{:<10s}: ({:6.4f}, {:6.4f}, {:6.4f}] to [{:<10s}:({:6.4f}, {:6.4f}, {:6.4f}]"
.format(
klist[i-1][3], klist[i-1][0], klist[i-1][1], klist[i-1][2],
klist[i][3], klist[i][0], klist[i][1], klist[i][2]))
dk_ = cal_kdistance(rg,
[klist[i-1][0], klist[i-1][1], klist[i-1][2]],
[klist[i][0], klist[i][1], klist[i][2]] )
ktotal += dk_
dklist.append(dk_)
ktotal_list.append(ktotal)
ktotal_namelist.append(klist[i][3])
# print(" dk={} ktotal={}".format(dk_, ktotal))
return dklist, ktotal_list, ktotal_namelist, ktotal
[ドキュメント]
def get_cal_klist(klist, nk):
"""
概要: バンド構造計算およびプロットのためのk点リストを生成します。
詳細説明:
`klist` と `nk` に基づき、k点パスに沿ってほぼ等間隔に配置された計算用k点ベクトル (`xkvec`) と、
プロット用の累積距離 (`xk`) を生成します。
また、k点パスの各セグメントに関する情報も辞書形式で返します。
:param klist: list k点パスのリスト。各要素は [kx, ky, kz, k点名] の形式。
:param nk: int プロットするk点の概数。
:returns xk: list プロット用のk点累積距離のリスト。
:returns xkvec: list 計算用のk点ベクトルリスト。各要素は [kx, ky, kz]。
:returns ktotallist: list k点パスの各高対称点における累積距離のリスト。
:returns ktotal_namelist: list k点パスの高対称点の名称リスト。
:returns res: dict k点計算に関する追加情報を含む辞書（"dklist", "nklist", "ktotal", "kstep"）。
"""
dklist, ktotallist, ktotal_namelist, ktotal = get_dklist(klist, nk)
kstep = ktotal / nk
nklist = []
xk = []
xkvec = []
ktotal_ = 0.0
for i in range(1, len(klist)):
nk_ = int(dklist[i-1] / kstep + 1.00001)
nklist.append(nk_)
if i == len(klist) - 1:
ndiv = nk_ - 1
else:
ndiv = nk_
kstepx_ = (klist[i][0] - klist[i-1][0]) / ndiv
kstepy_ = (klist[i][1] - klist[i-1][1]) / ndiv
kstepz_ = (klist[i][2] - klist[i-1][2]) / ndiv
print("k: {:<10s} - {:<10s}".format(klist[i-1][3], klist[i][3]), end = '')
print(": nk={:3d} ktotal={:6.4f} kstep=({:6.4f}, {:6.4f}, {:6.4f})"
.format(nklist[i-1], ktotal_, kstepx_, kstepy_, kstepz_))
dk_ = cal_kdistance(rg, [0.0, 0.0, 0.0], [kstepx_,kstepy_, kstepz_])
for j in range(nklist[i-1]):
kx = klist[i-1][0] + j * kstepx_
ky = klist[i-1][1] + j * kstepy_
kz = klist[i-1][2] + j * kstepz_
xk.append(ktotal_)
xkvec.append([kx, ky, kz])
ktotal_ += dk_
res = {"dklist": dklist, "nklist": nklist, "ktotal": ktotal, "kstep": kstep}
return xk, xkvec, ktotallist, ktotal_namelist, res
[ドキュメント]
def get_cal_Elist(xkvec, hrange, krange, lrange):
"""
概要: 指定されたk点リストと逆格子ベクトル範囲で自由電子エネルギーを計算します。
詳細説明:
各k点ベクトル `xkvec` に対して、指定された `hrange`, `krange`, `lrange` の
すべての逆格子ベクトルGを考慮した自由電子エネルギー `E(k+G)` を計算し、
リストとして返します。
:param xkvec: list 計算するk点ベクトル (kx, ky, kz) のリスト。
:param hrange: list H指数 [H_min, H_max]。
:param krange: list K指数 [K_min, K_max]。
:param lrange: list L指数 [L_min, L_max]。
:returns yE: list 各k点における複数のエネルギー値を含む入れ子のリスト。
"""
yE = []
for i in range(len(xkvec)):
kx = xkvec[i][0]
ky = xkvec[i][1]
kz = xkvec[i][2]
Elist = []
for ih in range(hrange[0], hrange[1]+1):
for ik in range(krange[0], krange[1]+1):
for il in range(lrange[0], lrange[1]+1):
E = cal_E([kx, ky, kz], [ih, ik, il])
Elist.append(E)
yE.append(Elist)
return yE
[ドキュメント]
def plot_band(axis, xk, yE, Erange, ktotallist, ktotal_namelist):
"""
概要: 自由電子バンド構造をプロットします。
詳細説明:
matplotlibのAxesオブジェクトに、計算されたバンドエネルギー `yE` をk点累積距離 `xk` に対してプロットします。
高対称点 (`ktotallist`, `ktotal_namelist`) の位置に縦線とラベルを追加し、エネルギー範囲を設定します。
:param axis: matplotlib.axes.Axes バンド構造をプロットするAxesオブジェクト。
:param xk: list プロットするk点の蓄積距離のリスト。
:param yE: list E(xk[i])。各k点におけるエネルギー値のリスト（入れ子になったリスト）。
:param Erange: list プロットするエネルギーの最小値と最大値 [E_min, E_max]。
:param ktotallist: list k点境界（高対称点）における、最初のk点からの距離の和のリスト。
:param ktotal_namelist: list k点境界における、k点の名称のリスト。
:returns: None
"""
# 表示範囲は決め打ち
axis.set_xlim([min(xk), max(xk)])
axis.set_ylim(Erange)
# バンド構造をプロット
axis.plot(xk, yE, linestyle = 'none',
marker = 'o', markerfacecolor = 'none', markeredgecolor = 'black',
markeredgewidth = 0.5, markersize = 2.0)
# Γ点、BZ境界の縦線を引く
for i in range(1, len(ktotallist)-1):
axis.plot([ktotallist[i], ktotallist[i]], Erange,
linestyle = '-', color = 'black', linewidth = 0.5)
# k軸の目盛りにk点の名称を表示する
# グラフ枠が一つであれば plt.xtics()で設定できる
# axisに対しては、.setpでattributeを直接書き換える必要があるらしい
plt.setp(axis, xticks = ktotallist, xticklabels = ktotal_namelist)
axis.set_xlabel("k", fontsize = fontsize)
axis.set_ylabel("E (eV)", fontsize = fontsize)
axis.tick_params(labelsize = fontsize)
[ドキュメント]
def main():
"""
概要: 自由電子バンド構造計算のメイン処理を実行します。
詳細説明:
結晶定数の表示、k点パスの処理、エネルギー計算、およびバンド構造のプロットを行います。
計算結果を標準出力に表示し、matplotlibによるプロットウィンドウを表示します。
ユーザーの入力があるまでウィンドウは閉じません。
"""
global a, rg
global nk, klist
print("")
print("Lattice parameters: ({}) A".format(a))
print("Effective mass: {} me".format(meff))
print("Reciprocal lattice metric [A^-2]:")
pprint(rg)
print("khl range: h=[{}, {}] k=[{}, {}] l=[{}, {}]"
.format(*hrange, *krange, *lrange))
print("Plot E range: {} - {} eV".format(*Erange))
xk, xkvec, ktotallist, ktotal_namelist, res = get_cal_klist(klist, nk)
print("")
print("k vectors")
print(" k_total: {} A^-1".format(res["ktotal"]))
print(" nk : ", nk)
print(" kstep : {}".format(res["kstep"]))
print(" dklist")
pprint(res["dklist"])
print(" ktotallist:")
pprint(ktotallist)
print(" nk_list:", res["nklist"])
# print("xk")
# pprint(xk)
yE = get_cal_Elist(xkvec, hrange, krange, lrange)
print("")
print("plot")
fig = plt.figure(figsize = figsize)
ax1 = fig.add_subplot(1, 1, 1)
plot_band(ax1, xk, yE, Erange, ktotallist, ktotal_namelist)
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
terminate()
if __name__ == "__main__":
main()