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
"""
1D band calculation by plain wave basis set
概要:
平面波基底セットを用いた1次元バンド計算を行うスクリプト。
詳細説明:
本スクリプトは、指定されたポテンシャルに対する1次元結晶の電子状態を、
平面波基底展開法によって計算します。
以下の3つのモードをサポートします。
- 'ft': ポテンシャルのフーリエ変換を表示します。
- 'band': バンド構造（エネルギーバンド図）を計算・表示します。
- 'wf': 特定のk点における波動関数と確率密度を計算・表示します。
コマンドライン引数により、計算モードや物理パラメータを変更可能です。
関連リンク:
:doc:`pw1d_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";
KE = hbar * hbar / 2.0 / me / e * 1.0e20 # coefficient of kinetic energy
#print("KE = hbar*hbar/(2m)/e*1.0e20 = {}", KE)
#========================
# global configuration
#========================
mode = 'ft' # ft|band|wf
#========================
# Crystal definition
#========================
# Si
a = 5.4064 # angstrom, lattice parameter
na = 64 # division for FFT, must be 2^n
#========================
# Potential
#========================
pottype = 'rect'
#pottype = 'gauss'
bwidth = 0.5 # A, barrier width
bpot = 10.0 # eV, barrier height
#========================
# Band
#========================
nG = 3 # # of G points (# of basis functions)
kmin = -0.5 # in pi/a
kmax = 0.5 # in pi/a
nk = 21
# プロットするエネルギー範囲
Erange = [0.0, 10.0] # eV
#========================
# Wave function
#========================
xwmin = 0.0 # A
xwmax = 3.0 * a # A
nxw = 101
kw = 0.0
iLevel = 0
#===================================
# figure configuration
#===================================
figsize = (6, 8)
fontsize = 12
legend_fontsize = 12
#==============================================
# fundamental functions
#==============================================
[ドキュメント]
def pfloat(str):
"""
概要: 文字列を浮動小数点数に変換する。
詳細説明:
文字列をfloat型に変換します。変換に失敗した場合はNoneを返しますが、
プログラムは終了しません。
Parameters:
:param str: str: 変換対象の文字列。
Returns:
:returns: float or None: 変換された浮動小数点数、または変換失敗時はNone。
"""
try:
return float(str)
except:
return None
[ドキュメント]
def pint(str):
"""
概要: 文字列を整数に変換する。
詳細説明:
文字列をint型に変換します。変換に失敗した場合はNoneを返しますが、
プログラムは終了しません。
Parameters:
:param str: str: 変換対象の文字列。
Returns:
:returns: int or None: 変換された整数、または変換失敗時はNone。
"""
try:
return int(str)
except:
return None
[ドキュメント]
def getarg(position, defval = None):
"""
概要: コマンドライン引数を取得する。
詳細説明:
sys.argvリストから指定された位置の引数を取得します。
指定された位置に引数が存在しない場合は、デフォルト値を返します。
Parameters:
:param position: int: 取得する引数のsys.argvリストにおけるインデックス。
:param defval: Any, optional: 引数が存在しない場合に返すデフォルト値。デフォルトはNone。
Returns:
:returns: str or Any: 取得された引数の文字列、またはデフォルト値。
"""
try:
return sys.argv[position]
except:
return defval
[ドキュメント]
def getfloatarg(position, defval = None):
"""
概要: コマンドライン引数を浮動小数点数として取得する。
詳細説明:
getarg()で取得したコマンドライン引数をpfloat()を用いて浮動小数点数に変換して返します。
Parameters:
:param position: int: 取得する引数のsys.argvリストにおけるインデックス。
:param defval: Any, optional: 引数が存在しない場合に返すデフォルト値。デフォルトはNone。
Returns:
:returns: float or None: 変換された浮動小数点数、またはNone。
"""
return pfloat(getarg(position, defval))
[ドキュメント]
def getintarg(position, defval = None):
"""
概要: コマンドライン引数を整数として取得する。
詳細説明:
getarg()で取得したコマンドライン引数をpint()を用いて整数に変換して返します。
Parameters:
:param position: int: 取得する引数のsys.argvリストにおけるインデックス。
:param defval: Any, optional: 引数が存在しない場合に返すデフォルト値。デフォルトはNone。
Returns:
:returns: int or None: 変換された整数、またはNone。
"""
return pint(getarg(position, defval))
[ドキュメント]
def usage():
"""
概要: プログラムの正しい使用方法を表示する。
詳細説明:
コマンドライン引数で異なるモード（ft, band, wf）を実行するための
正しい書式と例を標準出力に表示します。
"""
print("")
print("Usage: Variables in () are optional")
print(" python {}".format(sys.argv[0]))
print(" python {} (ft a na pottype bwidth bpot)".format(sys.argv[0]))
print(" python {} (band a na pottype bwidth bpot nG kmin kmax nk)".format(sys.argv[0]))
print(" python {} (wf a na pottype bwidth bpot nG kw iLevel xwmin xwmax nxw)".format(sys.argv[0]))
print(" pottype: rect|gauss")
print(" ex: python {} {} {} {} {} {} {}"
.format(sys.argv[0], 'ft', a, na, pottype, bwidth, bpot))
print(" ex: python {} {} {} {} {} {} {} {} {} {} {}"
.format(sys.argv[0], 'band', a, na, pottype, bwidth, bpot, nG, kmin, kmax, nk))
print(" ex: python {} {} {} {} {} {} {} {} {} {} {} {} {}"
.format(sys.argv[0], 'wf', a, na, pottype, bwidth, bpot, nG, kw, iLevel, xwmin, xwmax, nxw))
[ドキュメント]
def terminate(message = None):
"""
概要: メッセージを表示し、プログラムを終了する。
詳細説明:
オプションのメッセージを表示した後、usage()関数を呼び出して使用方法を表示し、
プログラムを終了します。
Parameters:
:param message: str, optional: 終了時に表示するメッセージ。デフォルトはNone。
"""
print("")
if message is not None:
print("")
print(message)
print("")
usage()
print("")
exit()
#==============================================
# update default values by startup arguments
#==============================================
argv = sys.argv
#if len(argv) == 1:
# terminate()
mode = getarg (1, mode)
a = getfloatarg(2, a)
na = getintarg (3, na)
pottype = getarg (4, pottype)
bwidth = getfloatarg(5, bwidth)
bpot = getfloatarg(6, bpot)
if mode == 'band':
nG = getintarg ( 7, nG)
kmin = getfloatarg( 8, kmin)
kmax = getfloatarg( 9, kmax)
nk = getintarg (10, nk)
elif mode == 'wf':
nG = getintarg ( 7, nG)
kw = getfloatarg( 8, kw)
iLevel = getintarg ( 9, iLevel)
xwmin = getfloatarg(10, xwmin)
xwmax = getfloatarg(11, xwmax)
nxw = getintarg (12, nxw)
[ドキュメント]
def reduce(x, x0):
"""
概要: 空間座標xを結晶格子定数x0で周期化し、基本セル内の座標に変換する。
詳細説明:
与えられた座標xを、周期x0の範囲（例: 0からx0）に収まるように変換します。
これにより、周期的なポテンシャルや波動関数の計算で、基本セル内での
座標を扱うことができます。
Parameters:
:param x: float: 元の空間座標。
:param x0: float: 周期的な長さ（例: 格子定数a）。
Returns:
:returns: float: 基本セル内に正規化された空間座標。
"""
n = int(x / x0)
if x < 0.0:
n += 1
return x - x0 * n
[ドキュメント]
def pot(x):
"""
概要: 1次元の周期的なポテンシャル値を計算する。
詳細説明:
グローバル変数`pottype`('rect'または'gauss')、`bwidth`、`bpot`、`a`
に基づいて、指定された空間座標xにおけるポテンシャル値を計算します。
ポテンシャルは格子定数`a`で周期化されます。
Parameters:
:param x: float: 空間座標。
Returns:
:returns: float: 空間座標xにおけるポテンシャル値 (eV)。
"""
global a
global pottype, bwidth, bpot
xred = reduce(x, a)
if pottype == 'rect':
if 0.0 <= xred <= bwidth:
return bpot
return 0.0
if pottype == 'gauss':
xx = (xred - 0.5 * a) / bwidth
return bpot * exp(-xx * xx)
[ドキュメント]
def build_potential(xmin, xstep, n):
"""
概要: 指定された範囲とステップでポテンシャルプロファイルを作成する。
詳細説明:
xminから開始し、xstep間隔でn個の点におけるポテンシャル値を計算し、
空間座標とポテンシャル値の配列を返します。
Parameters:
:param xmin: float: ポテンシャルプロファイルの開始座標。
:param xstep: float: 座標のステップサイズ。
:param n: int: 計算する点の数。
Returns:
:returns: tuple[numpy.ndarray, numpy.ndarray]:
- xpot (numpy.ndarray): 各点の空間座標の配列。
- ypot (numpy.ndarray): 各点におけるポテンシャル値の配列。
"""
xpot = np.empty(n)
ypot = np.empty(n)
for i in range(n):
xx = xmin + i * xstep
xpot[i] = xx
ypot[i] = pot(xx)
return xpot, ypot
[ドキュメント]
def find_iG(dij, Glist):
"""
概要: 逆格子点リストの中から指定された逆格子内部座標に対応するインデックスを検索する。
詳細説明:
Glistは、平面波基底として使用される逆格子点の内部座標（整数）のリストです。
dijがこのリスト内に存在する場合、そのインデックスを返します。
Parameters:
:param dij: int: 検索対象の逆格子内部座標。
:param Glist: numpy.ndarray: 逆格子点内部座標のリスト。
Returns:
:returns: int or None: dijに対応するGlist内のインデックス、または見つからない場合はNone。
"""
for iG in range(len(Vft)): # Vftはグローバル変数ではなくVftlistであるべきだが、既存コードを尊重
if Glist[iG] == dij:
return iG
return None
[ドキュメント]
def find_Vft(dij, Glist, Vft):
"""
概要: 指定された逆格子内部座標に対応するフーリエ変換されたポテンシャル成分を返す。
詳細説明:
逆格子点Glistと対応するフーリエ変換ポテンシャルVftのリストから、
指定された逆格子内部座標dijに対応するVft値を検索して返します。
見つからない場合は0.0を返します。
Parameters:
:param dij: int: 検索対象の逆格子内部座標。
:param Glist: numpy.ndarray: 逆格子点内部座標のリスト。
:param Vft: numpy.ndarray: 各G点に対応するフーリエ変換されたポテンシャル成分。
Returns:
:returns: complex: dijに対応するフーリエ変換されたポテンシャル成分、または0.0。
"""
for iG in range(len(Vft)):
if Glist[iG] == dij:
return Vft[iG]
return 0.0
[ドキュメント]
def cal_fft(na, a, ypot):
"""
概要: 1次元ポテンシャルのフーリエ変換を計算し、物理的なG点座標に変換する。
詳細説明:
与えられたポテンシャル配列ypotに対して、numpy.fft.fftを用いてフーリエ変換を行います。
FFTの結果は、負のG成分が配列の後半にあるため、物理的なG点座標（負から正）に対応するように
再配置・変換します。
Parameters:
:param na: int: FFTに用いる点の数（通常は2のN乗）。
:param a: float: 格子定数 (Å)。
:param ypot: numpy.ndarray: 空間座標におけるポテンシャル値の配列。
Returns:
:returns: tuple[list, numpy.ndarray, numpy.ndarray, numpy.ndarray, list, int, float]:
- xft0 (list): raw FFT結果のインデックス。
- yft0 (numpy.ndarray): raw FFT結果の複素数配列。
- xft (numpy.ndarray): 物理的なG点座標（単位1/Å）の配列。
- yft (numpy.ndarray): 物理的なG点座標に対応するフーリエ変換ポテンシャルの複素数配列。
- iGlist (list): 各FFT成分に対応する逆格子点内部座標 (iG) のリスト。
- nahalf (int): naの半分。
- xftstep (float): G点座標のステップサイズ (1/Å)。
"""
xftstep = 1.0 / a
astep = a / na
nahalf = int(na/2)
xftmax = nahalf * xftstep
xft0 =range(na)
yft0 = np.fft.fft(ypot) * astep / a
# note: The fft result has periodicity
# The second half at i >= nahalf data should be shifted to negative x
xft = np.arange(-xftmax, xftmax, xftstep)
yft = np.empty(na, dtype = complex)
for i in range(nahalf):
yft[i] = yft0[nahalf + i]
for i in range(nahalf):
yft[i+nahalf] = yft0[i]
iGlist = []
for i in range(na):
if i <= nahalf:
iG = i
else:
iG = -(na - i)
iGlist.append(iG)
return xft0, yft0, xft, yft, iGlist, nahalf, xftstep
[ドキュメント]
def free_KE(k, iG):
"""
概要: 自由電子の運動エネルギーを計算する。
詳細説明:
電子の運動エネルギーを、波数kと逆格子点内部座標iGの合計を用いて計算します。
単位はeVです。
Parameters:
:param k: float: 電子の波数（k点、単位はπ/a）。
:param iG: int: 逆格子点の内部座標（整数）。
Returns:
:returns: float: 自由電子の運動エネルギー (eV)。
"""
kK = KE * pow(2.0*pi/a, 2.0) # coefficient of kinetic energy in eV and angstrom
return kK * (k + iG) * (k + iG)
[ドキュメント]
def solve_pw(Glist, Vftlist, nGmax, nG, k, IsPrint = 0):
"""
概要: 平面波基底を用いた固有方程式を解き、固有エネルギーと波動関数の係数を計算する。
詳細説明:
指定されたG点リストとフーリエ変換されたポテンシャルを用いて、
平面波基底によるハミルトニアン行列（Fock行列）を構築します。
この行列を対角化することで、固有エネルギー（バンド構造のエネルギー準位）と
それに対応する波動関数の平面波係数を計算します。
波動関数は規格化されていません。
Parameters:
:param Glist: numpy.ndarray: 逆格子点内部座標の配列。
:param Vftlist: numpy.ndarray: 各G点に対応するフーリエ変換されたポテンシャル成分の配列。
:param nGmax: int: 抽出されたG点の最大絶対値。
:param nG: int: 基底関数の数。
:param k: float: 電子の波数（k点、単位はπ/a）。
:param IsPrint: int, optional: デバッグ情報を出力するレベル。0:なし, 1:一部, 2:詳細。デフォルトは0。
Returns:
:returns: tuple[numpy.ndarray, numpy.ndarray]:
- ei (numpy.ndarray): 計算された固有エネルギーの配列 (eV)。
- ci (numpy.ndarray): 各固有エネルギーに対応する平面波係数の行列。
ci[i][j]はi番目のG点におけるj番目のバンドの係数。
"""
kK = KE * pow(2.0*pi/a, 2.0) # coefficient of kinetic energy in eV and angstrom
# build Fock matrix
Hij = np.empty([nG, nG], dtype = complex)
for i in range(nG):
iG = Glist[i]
# Diagonal elements: Kinetic Energy + V(G=0)
# find_Vft(0, Glist, Vftlist) はG=0の成分を取得する。
# VftlistはGlistに対応するフーリエ係数のリストなので、Glist[0]が0であると仮定できる。
# 実際にはextract_basisでGlist[0]=0, Vftlist[0]=yft0[0]となっている。
Vij = Vftlist[0] # find_Vft(0, Glist, Vftlist)
Hij[i][i] = free_KE(k, iG) + Vij
for j in range(i+1, nG):
jG = Glist[j]
dij = iG - jG # 逆格子ベクトル間の差
Vij = find_Vft(dij, Glist, Vftlist) # ポテンシャルのフーリエ成分 V(G_i - G_j)
if IsPrint >= 1:
print(" iG,jG,dij=", iG, jG, dij, " nG(base)=[{}, {}]".format(-nGmax, nGmax), " Vij=", Vij)
Hij[i][j] = Vij
Hij[j][i] = Hij[i][j].conjugate()
# print(" Hij=\n", Hij)
ei, ci = LA.eig(Hij)
if IsPrint >= 1:
print(" ei=", ei)
if IsPrint >= 2:
print(" ci=", ci)
return ei, ci
[ドキュメント]
def cal_wf(xwmin, xwmax, nxw, kw, ci, Glist):
"""
概要: 平面波基底の係数から波動関数と確率密度を計算する。
詳細説明:
指定されたk点、平面波基底の係数ci、およびG点リストを用いて、
波動関数（複素数）と、その絶対値の2乗である確率密度（実数）を、
指定された空間範囲で計算します。波動関数は規格化されていません。
Parameters:
:param xwmin: float: 波動関数の計算を開始する空間座標 (Å)。
:param xwmax: float: 波動関数の計算を終了する空間座標 (Å)。
:param nxw: int: 波動関数を計算する点の数。
:param kw: float: 電子の波数（k点、単位はπ/a）。
:param ci: list: 選択されたバンドの平面波係数のリスト。ci[iG] はGlist[iG]に対応する係数。
:param Glist: numpy.ndarray: 逆格子点内部座標の配列。
Returns:
:returns: tuple[numpy.ndarray, numpy.ndarray, float]:
- xwf (numpy.ndarray): 波動関数が計算された空間座標の配列。
- ywf (numpy.ndarray): 各空間点における波動関数（複素数）の配列。
- charge (float): 最後の点の確率密度。
※この戻り値は関数内で累積されていないため、正確には最後の点の確率密度を返す。
実際にはywf2[i]で各点の確率密度が計算されているが、最後の値しか返されない。
コードの変更不可ルールによりこのまま。
"""
nG = len(ci)
xwstep = (xwmax - xwmin) / (nxw - 1)
xwf = np.empty(nxw, dtype = complex) # 実際には実数
ywf = np.empty(nxw, dtype = complex)
ywf2 = np.empty(nxw, dtype = float)
for i in range(nxw):
x = xwmin + i * xwstep
f = 0.0 + 0.0j # 複素数として初期化
for iG in range(nG):
G = Glist[iG]
f += ci[iG] * exp(1.0j * (kw+G) * pi2/a * x)
# print("f({})={}".format(x, f))
xwf[i] = x
ywf[i] = f
charge = f * f.conjugate()
ywf2[i] = charge.real
return xwf, ywf, charge
[ドキュメント]
def wf():
"""
概要: 波動関数計算モードのメイン処理を実行する。
詳細説明:
グローバル変数で定義された結晶、ポテンシャル、k点、エネルギー準位、空間範囲に基づき、
波動関数とその確率密度を計算し、プロットします。
ポテンシャルのフーリエ変換、固有方程式の解法、波動関数の計算と可視化が含まれます。
"""
global mode
global a, na
global kw, iG
global xwmin, xwmax, nxw
xwstep = (xwmax - xwmin) / (nxw - 1)
print("")
print("=== Input parameterss ===")
print("mode:", mode)
print("a=", a, "A")
print(" na=", na)
print("potential: {} w={} A h={} eV".format(pottype, bwidth, bpot))
print("Wave function to be plotted: k = {} iLevel = {}".format(kw, iLevel))
print("x range: {} - {} at {} step, {} points".format(xwmin, xwmax, xwstep, nxw))
print("potential: {} w={} A h={} eV".format(pottype, bwidth, bpot))
astep = a / na
xpot, ypot = build_potential(0.0, astep, na)
xplot, yplot = build_potential(xwmin, xwstep, nxw)
xft0, yft0, xft, yft, iGlist, nahalf, xftstep = cal_fft(na, a, ypot)
print("")
print("=== FT result ===")
print("{:4} {:4} {}".format("i", "iG", "ci"))
for i in range(na):
print("{:4d} {:4d} {}".format(i, iGlist[i], yft0[i]))
# Extract G basis set for given nG
Glist, Vftlist, nGmax = extract_basis(yft0, nG)
print("")
print("=== G basis extracted ===")
print("nG =", nG)
print("{:4} {}".format("iG", "Vft"))
for i in range(nG):
print("{:4d} {}".format(Glist[i], Vftlist[i]))
print("")
print("=== Solve eigen equations ===")
k = kw
print("at k = {}".format(k))
yE, ci = solve_pw(Glist, Vftlist, nGmax, nG, k, IsPrint = 1)
print(" ei=", yE)
print(" ci=", ci)
print("")
print("=== Calculate wave function ===")
print("Energy levels:")
for i in range(len(yE)):
print(" {} {:12.6g} eV".format(i, yE[i].real))
print("")
print("at k = {}".format(k))
print(" selected for {}-th energy level:".format(iLevel))
print(" E = {:12.6g} eV".format(yE[iLevel].real))
print(" Wave function coefficinets")
# ciは (nG, nG) の行列で、ci[i][j] は i番目の基底 Glist[i] に対する j番目の固有状態の係数
ciLevel = [ci[i][iLevel] for i in range(nG)] # iLevel番目の固有状態の係数を抽出
for i in range(nG):
print(" at iG={:3d}".format(Glist[i]), " ci={:8.4f}".format(ci[i][iLevel]))
xwf, ywf, charge_at_last_point = cal_wf(xwmin, xwmax, nxw, kw, ciLevel, Glist) # chargeは最後の点の値のみ返されるため注意
charge = [(f0 * f0.conjugate()).real for f0 in ywf] # 全ての点の確率密度を再計算
print("c")
fig = plt.figure(figsize = (8, 8))
ax2 = fig.add_subplot(1, 1, 1)
# ax2 = fig.add_subplot(2, 1, 2)
# ax2をax1に関連させる
ax1 = ax2.twinx()
ax1.set_xlim([xwmin, xwmax])
ax1.plot(xplot, yplot, linewidth = 0.5, label = '$V$($x$)')
ax1.plot(ax1.get_xlim(), [0.0, 0.0], color = 'r', linestyle = 'dashed', linewidth = 0.5)
ax2.set_xlim([xwmin, xwmax])
ax2.plot(xplot, ywf.real, linewidth = 1.5, label = r"$\Psi$(real)")
ax2.plot(xplot, ywf.imag, linewidth = 1.5, label = r"$\Psi$(imaginary)")
ax2.plot(xplot, charge, linewidth = 0.5, label = r"|$\Psi$|$^2$")
ax2.plot(ax1.get_xlim(), [0.0, 0.0], color = 'r', linestyle = 'dashed', linewidth = 0.5)
# ax1.set_xlabel(r"$x$ ($\AA$)", fontsize = fontsize)
ax1.set_ylabel(r"$V$($x$)", fontsize = fontsize)
ax2.set_xlabel(r"$x$ ($\AA$)", fontsize = fontsize)
ax2.set_ylabel(r"$\Psi$", fontsize = fontsize)
# 凡例をまとめて出力する
handler1, label1 = ax1.get_legend_handles_labels()
handler2, label2 = ax2.get_legend_handles_labels()
ax2.legend(handler1 + handler2, label1 + label2, loc = 2, borderaxespad = 0.0, fontsize = legend_fontsize)
# ax2.legend(fontsize = legend_fontsize)
ax1.tick_params(labelsize = fontsize)
ax2.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
terminate()
[ドキュメント]
def band():
"""
概要: バンド構造計算モードのメイン処理を実行する。
詳細説明:
グローバル変数で定義された結晶、ポテンシャル、G点基底の数、k点範囲に基づき、
1次元のエネルギーバンド構造を計算し、プロットします。
自由電子のエネルギーバンドも比較のために表示されます。
"""
global mode
global a, na
global nG
global kmin, kmax, nk
global pottype, bwidth, bpot
kstep = (kmax - kmin) / (nk - 1)
print("")
print("=== Input parameterss ===")
print("mode:", mode)
print("a=", a, "A")
print(" na=", na)
print("potential: {} w={} A h={} eV".format(pottype, bwidth, bpot))
print("Basis function: nG=", nG)
if nG >= 4:
nG = int(nG / 2) * 2 + 1 # convert to an odd number for nG >= 4
elif nG <= 0 or 64 < nG:
terminate("Error: nG must be between 1 and 64 [nG={}]".format(nG))
print("k range: {} - {} at {} step, {} points".format(kmin, kmax, kstep, nk))
print("potential: {} w={} A h={} eV".format(pottype, bwidth, bpot))
astep = a / na
xpot, ypot = build_potential(0.0, astep, na)
xft0, yft0, xft, yft, iGlist, nahalf, xftstep = cal_fft(na, a, ypot)
print("")
print("=== FT result ===")
print("{:4} {:4} {}".format("i", "iG", "ci"))
for i in range(na):
print("{:4d} {:4d} {}".format(i, iGlist[i], yft0[i]))
# Extract G basis set for given nG
Glist, Vftlist, nGmax = extract_basis(yft0, nG)
print("")
print("=== FTed potential ===")
print("nG =", nG)
print("{:4} {}".format("iG", "Vft"))
for i in range(nG):
print("{:4d} {}".format(Glist[i], Vftlist[i]))
print("")
print("=== Solve eigen equations ===")
xk = [kmin + i * kstep for i in range(nk)]
yE = np.empty([nG, nk])
yEfree = np.empty([nG, nk])
for ik in range(nk):
k = kmin + ik * kstep
print("at k = {}".format(k))
ei,ci = solve_pw(Glist, Vftlist, nGmax, nG, k, IsPrint = 1)
for i in range(nG):
yE[i][ik] = ei[i]
# yEfreeはGlistに対応する自由電子バンドを計算
# ただし、eiがソートされている場合、iGはGlist[i]とは一致しない可能性があるため注意
# ここでは各G点に対応する自由電子のエネルギーを計算し、nG個のバンドとしてプロットしている。
# プロット時にはyEもyEfreeも各G点に対してプロットされるため、見かけ上問題ない。
# (Glist[0] = 0, Glist[1] = 1, Glist[2] = -1 のような順序で格納されている)
iG = Glist[i]
print("k=", k, iG, k+iG)
yEfree[i][ik] = free_KE(k, iG)
fig = plt.figure(figsize = figsize)
ax1 = fig.add_subplot(1, 1, 1)
ax1.set_xlim([-0.5, 0.5])
ax1.set_ylim(Erange)
# yEfree はnG個の自由電子バンドを描画 (Glistの各要素に対する自由電子)
for iG_idx in range(nG): # Glistのインデックス
ax1.plot(xk, yE[iG_idx], linestyle = '', marker = 'o', markersize = 2.0,
markerfacecolor = 'none', markeredgecolor = 'black', markeredgewidth = 0.5) #, label = 'with V(x)')
ax1.plot(xk, yEfree[iG_idx], linestyle = '-', linewidth = 0.3, color = 'red')# , label = 'free e')
ax1.set_xlabel(r"k ($\pi$/a)", fontsize = fontsize)
ax1.set_ylabel("E (eV)", fontsize = fontsize)
# 現状のコードでは凡例は表示されない。label = 'with V(x)'と 'free e' がコメントアウトされているため。
ax1.legend(fontsize = legend_fontsize)
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
terminate()
[ドキュメント]
def ft():
"""
概要: ポテンシャルのフーリエ変換表示モードのメイン処理を実行する。
詳細説明:
グローバル変数で定義された結晶とポテンシャルに基づいて、
空間座標におけるポテンシャル、およびそのフーリエ変換結果をプロットします。
フーリエ変換の生データと、物理的なG点座標に変換されたデータが両方表示されます。
"""
global mode
global a, na
global pottype, bwidth, bpot
astep = a / na
print("")
print("=== Input parameterss ===")
print("mode:", mode)
print("a=", a, "A")
print(" na=", na)
print(" astep=", astep)
print("potential: {} w={} A h={} eV".format(pottype, bwidth, bpot))
print("")
xpot, ypot = build_potential(0.0, astep, na)
xplot, yplot = build_potential(0.0, astep, 3 * na)
xft0, yft0, xft, yft, iG, nahalf, xftstep = cal_fft(na, a, ypot)
fig = plt.figure(figsize = (8, 8))
ax1 = fig.add_subplot(2, 2, 1)
ax2 = fig.add_subplot(2, 2, 3)
ax3 = fig.add_subplot(2, 2, 4)
ax1.plot(xplot, yplot, linewidth = 0.5)
ax2.plot(xft0, yft0.real, label = "real", linewidth = 0.5)
ax2.plot(xft0, yft0.imag, label = "imaginary", linewidth = 0.5)
ax3.plot(xft, yft.real, marker = 'o', markersize = 1.0, linewidth = 0.5, label = "real")
ax3.plot(xft, yft.imag, marker = 'o', markersize = 1.0, linewidth = 0.5, label = "imaginary")
ax3.plot(xft, abs(yft), marker = 'o', markersize = 1.0, linewidth = 0.5, label = "absolute")
ax1.set_xlabel(r"x ($\AA$)", fontsize = fontsize)
ax1.set_ylabel("pot (x)", fontsize = fontsize)
ax2.set_xlabel("i", fontsize = fontsize)
ax2.set_ylabel("FTed pot, raw x", fontsize = fontsize)
ax3.set_xlabel("1/x, normalized (1/A)", fontsize = fontsize)
ax3.set_ylabel("FTed pot")
ax3.set_xlim([-10.0 * xftstep, 10.0 * xftstep])
# ax1.legend(fontsize = legend_fontsize)
ax2.legend(fontsize = legend_fontsize)
ax3.legend(fontsize = legend_fontsize)
ax1.tick_params(labelsize = fontsize)
ax2.tick_params(labelsize = fontsize)
ax3.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
terminate()
[ドキュメント]
def main():
"""
概要: プログラムのメインエントリーポイント。
詳細説明:
グローバル変数`mode`の値に基づいて、適切な計算モードの関数を呼び出します。
サポートされるモードは 'ft' (フーリエ変換表示), 'band' (バンド構造計算),
'wf' (波動関数計算) です。認識できないモードの場合はプログラムを終了します。
"""
global mode
if mode == 'ft':
ft()
elif mode == 'band':
band()
elif mode == 'wf':
wf()
else:
terminate("Error: Invalid mode [{}]".format(mode))
if __name__ == "__main__":
main()