"""
このスクリプトは、一次元の周期ポテンシャルにおける電子のバンド構造や波動関数を平面波基底を用いて計算します。
概要:
フーリエ変換されたポテンシャルを利用してハミルトニアン行列を構築し、固有方程式を解くことで
エネルギーバンドや波動関数係数を求めます。
詳細説明:
1. 周期ポテンシャル（長方形またはガウス型）を空間領域で構築します。
2. FFTを用いてポテンシャルを逆格子空間へ変換します。
3. 指定された数の平面波基底（G点）を用いて、各波数 k におけるハミルトニアン行列を構築します。
4. 対角化により固有エネルギー（バンド）と固有ベクトル（波動関数係数）を算出します。
5. 得られた係数を用いて、実空間での波動関数および電荷密度を合成・可視化します。
関連リンク: :doc:`pw1d_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
#===================================
# 物理定数
#===================================
PI = 3.14159265358979323846
PI2 = 2.0 * PI
H_PLANCK = 6.6260755e-34 # Js
HBAR = 1.05459e-34 # Js
C_LIGHT = 2.99792458e8 # m/s
E_CHARGE = 1.60218e-19 # C
EPS0 = 8.854418782e-12 # C^2N^-1m^-2
KB = 1.380658e-23 # JK^-1
ME = 9.1093897e-31 # kg
A0 = 5.29177e-11 # m
# 運動エネルギー係数 (hbar^2 / 2m) [eV * A^2]
KE_COEFF = (HBAR**2) / (2.0 * ME * E_CHARGE) * 1.0e20
#==============================================
# ヘルパー関数群
#==============================================
[ドキュメント]
def pfloat(s):
"""文字列を浮動小数点数に安全に変換します。"""
try: return float(s)
except: return None
[ドキュメント]
def pint(s):
"""文字列を整数に安全に変換します。"""
try: return int(s)
except: return None
[ドキュメント]
def getarg(position, defval=None):
"""コマンドライン引数を安全に取得します。"""
try: return sys.argv[position]
except: return defval
[ドキュメント]
def reduce_range(x, x0):
"""実数値 x を周期 x0 の範囲 [0, x0) に収めます。"""
return x % x0
[ドキュメント]
def pot_func(x, pottype, a, bwidth, bpot):
"""
一次元周期ポテンシャルを計算します。
:param x: 位置 [A]
:param pottype: 'rect' または 'gauss'
:param a: 格子定数 [A]
"""
xred = reduce_range(x, a)
if pottype == 'rect':
return bpot if 0.0 <= xred <= bwidth else 0.0
if pottype == 'gauss':
xx = (xred - 0.5 * a) / bwidth
return bpot * exp(-xx * xx)
return 0.0
[ドキュメント]
def build_potential(xmin, xstep, n, pottype, a, bwidth, bpot):
"""指定された範囲でポテンシャルプロファイルを構築します。"""
xpot = np.linspace(xmin, xmin + xstep * (n-1), n)
ypot = np.array([pot_func(xx, pottype, a, bwidth, bpot) for xx in xpot])
return xpot, ypot
[ドキュメント]
def cal_fft(na, a, ypot):
"""ポテンシャルのフーリエ変換を計算し、逆格子空間の形式に整理します。"""
astep = a / na
yft0 = np.fft.fft(ypot) * astep / a
xftstep = 1.0 / a
nahalf = int(na/2)
xftmax = nahalf * xftstep
xft = np.arange(-xftmax, xftmax, xftstep)
yft = np.concatenate([yft0[nahalf:], yft0[:nahalf]])
iGlist = []
for i in range(na):
iG = i if i <= nahalf else -(na - i)
iGlist.append(iG)
return yft0, xft, yft, iGlist, xftstep
[ドキュメント]
def find_Vft(dij, iGlist, yft0):
"""指定された逆格子座標 dij に対応するフーリエ成分を返します。"""
# FFTの結果のインデックスに直接アクセス
idx = dij if dij >= 0 else len(yft0) + dij
if 0 <= idx < len(yft0):
return yft0[idx]
return 0.0
[ドキュメント]
def solve_pw(Glist, yft0, nG, k, a):
"""平面波基底を用いたハミルトニアンの対角化。"""
k_factor = KE_COEFF * (2.0 * PI / a)**2
Hij = np.zeros((nG, nG), dtype=complex)
for i in range(nG):
iG = Glist[i]
# 対角項: (k+G)^2 + V(0)
Hij[i, i] = k_factor * (k + iG)**2 + yft0[0]
for j in range(i + 1, nG):
jG = Glist[j]
# 非対角項: V(Gi-Gj)
Vij = find_Vft(iG - jG, None, yft0)
Hij[i, j] = Vij
Hij[j, i] = Vij.conjugate()
ei, ci = LA.eig(Hij)
# エネルギーの低い順にソート
idx = ei.argsort()
return ei[idx], ci[:, idx]
#==============================================
# モード別実行関数
#==============================================
[ドキュメント]
def run_band(config):
"""バンド構造の計算とプロット"""
a, na, nG = config['a'], config['na'], config['nG']
pottype, bwidth, bpot = config['pottype'], config['bwidth'], config['bpot']
kmin, kmax, nk = config['kmin'], config['kmax'], config['nk']
xpot, ypot = build_potential(0.0, a/na, na, pottype, a, bwidth, bpot)
yft0, xft, yft, iGlist, _ = cal_fft(na, a, ypot)
# G基底の抽出 (0, 1, -1, 2, -2...)
nG_actual = min(nG, na)
Glist = [i if i <= nG_actual//2 else -(nG_actual - i) for i in range(nG_actual)]
Glist.sort(key=abs)
Glist = Glist[:nG]
xk = np.linspace(kmin, kmax, nk)
yE = np.zeros((nG, nk))
yEfree = np.zeros((nG, nk))
k_factor = KE_COEFF * (2.0 * PI / a)**2
for ik, k in enumerate(xk):
ei, _ = solve_pw(Glist, yft0, nG, k, a)
for i in range(nG):
yE[i, ik] = ei[i].real
yEfree[i, ik] = k_factor * (k + Glist[i])**2
plt.figure(figsize=(6, 8))
for i in range(nG):
plt.plot(xk, yE[i], 'o', markersize=2, color='black', mfc='none')
plt.plot(xk, yEfree[i], '-', linewidth=0.3, color='red', alpha=0.5)
plt.ylim(0, 15)
plt.xlabel("k (pi/a)")
plt.ylabel("E (eV)")
plt.title(f"1D Band Structure ({pottype} potential)")
plt.show()
[ドキュメント]
def main():
# デフォルト設定
config = {
'mode': 'band', 'a': 5.4064, 'na': 64,
'pottype': 'rect', 'bwidth': 0.5, 'bpot': 10.0,
'nG': 11, 'kmin': -1.0, 'kmax': 1.0, 'nk': 41
}
# 引数による更新 (簡易実装)
if len(sys.argv) > 1:
config['mode'] = sys.argv[1]
if config['mode'] == 'band':
run_band(config)
else:
print("Selected mode is under construction or not implemented in this simplified version.")
if __name__ == "__main__":
main()