"""
1Dタイトバインディング法による量子井戸ポテンシャルを持つバンド計算スクリプト。
このスクリプトは、一次元タイトバインディング法を用いて、量子井戸ポテンシャル下での電子のバンド構造や波動関数を計算し、可視化します。
原子配列やポテンシャル形状、計算モード（basis, band, wf）などのパラメータを設定ファイル的に変更して実行できます。
:doc:`tb1d_usage`
"""
import sys
import numpy as np
from numpy import sqrt, exp, sin, cos, tan, sinh, cosh, tanh, 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";
KE = hbar * hbar / 2.0 / me / e * 1.0e20 # coefficient of kinetic energy
#========================
# global parameters
#========================
parity_str = ["even", "odd"]
#========================
# global configuration
#========================
mode = 'basis' # basis|band|wf
#========================
# Potential
#========================
Vvac = 10.0 # eV
#========================
# Crystal definition
#========================
#crystalname = 'Sim'
#crystalname = 'AB'
crystalname = 'Si'
# Si
a = 5.4064 # angstrom, lattice parameter
# tijを計算する近接サイトの検索範囲
rmax_neighbor = 1.1 * a
atoms = [
{"name": "A", "wellwidth": a - 0.5, "Vbase": 0.0, "meff": 1.0}
]
crystal = [
{"name": 'A1', "iatom": 0, "x0": (0.5 + a) / 2.0}
]
if crystalname == 'AB':
a = 5.4064 # angstrom, lattice parameter
rmax_neighbor = 1.1 * a
atoms = [
{"name": "A", "wellwidth": 2.5, "Vbase": 0.0, "meff": 5.0}
,{"name": "B", "wellwidth": a - 3.5, "Vbase": 0.3, "meff": 5.0}
]
crystal = [
{"name": 'A1', "iatom": 0, "x0": (0.5 + 3.0)/2.0}
,{"name": 'B1', "iatom": 1, "x0": (3.5 + a)/2.0}
]
if crystalname == 'Sim':
a = 17.5
rmax_neighbor = 4.0
atoms = [
{"name": "A", "wellwidth": 3.0, "Vbase": 0.0, "meff": 1.0}
]
crystal = [
{"name": 'A1', "iatom": 0, "x0": 2.0}
,{"name": 'A2', "iatom": 0, "x0": 5.5}
,{"name": 'A3', "iatom": 0, "x0": 9.0}
,{"name": 'A4', "iatom": 0, "x0": 12.5}
,{"name": 'A5', "iatom": 0, "x0": 16.0}
]
#==============================================
# X range to plot potential and basis
#==============================================
xmin = -0.1 * a # A
xmax = 1.1 * a # A
nx = 301
#=====================================
# 解を走査するグラフ表示
#=====================================
# 解を走査するエネルギー範囲
Emin = 0.0
Emax = Vvac
# 解を走査するエネルギー点数
nE = 101
# Secant法パラメータ
eps = 1.0e-5
nmaxiter = 100
#========================
# Wave function plot
#========================
wmargin = 0.5
nwf = 101
wfplotmax = 1.0 # superiposed on energy level diagram, in eV
#========================
# <E> integration
#========================
xminEinteg = -10.0
xmaxEinteg = a + 11.0
nxEinteg = 501
#========================
# Band
#========================
# tijを計算する近接サイトを検索する格子周期範囲
Nxmax_neighbor = 2
# tijを計算する際の波動関数の計算範囲
tijxmin = -1.1 * a
tijxmax = 1.1 * a
ntijx = 501
# band plot
kmin = -0.5 # in pi/a
kmax = 0.5 # in pi/a
nk = 21
# プロットするエネルギー範囲
Erange = [0.0, 10.0] # eV
#===================================
# figure configuration
#===================================
figsize = (6, 8)
fontsize = 12
legend_fontsize = 12
#==============================================
# fundamental functions
#==============================================
[ドキュメント]
def pfloat(str):
"""
文字列を浮動小数点数に安全に変換します。
変換できない場合はNoneを返します。
:param str: 変換する文字列。
:returns: 変換された浮動小数点数、またはNone。
"""
try:
return float(str)
except:
return None
[ドキュメント]
def pint(str):
"""
文字列を整数に安全に変換します。
変換できない場合はNoneを返します。
:param str: 変換する文字列。
:returns: 変換された整数、またはNone。
"""
try:
return int(str)
except:
return None
[ドキュメント]
def getarg(position, defval = None):
"""
コマンドライン引数を安全に取得します。
指定された位置に引数がない場合はデフォルト値を返します。
:param position: 取得する引数の位置（0から始まる）。
:param defval: 引数が見つからない場合に返すデフォルト値。
:returns: 指定された位置の引数、またはデフォルト値。
"""
try:
return sys.argv[position]
except:
return defval
[ドキュメント]
def getfloatarg(position, defval = None):
"""
コマンドライン引数を浮動小数点数として安全に取得します。
取得できない場合や変換できない場合はデフォルト値を返します。
:param position: 取得する引数の位置（0から始まる）。
:param defval: 引数が見つからない場合や変換できない場合に返すデフォルト値。
:returns: 変換された浮動小数点数、またはデフォルト値。
"""
return pfloat(getarg(position, defval))
[ドキュメント]
def getintarg(position, defval = None):
"""
コマンドライン引数を整数として安全に取得します。
取得できない場合や変換できない場合はデフォルト値を返します。
:param position: 取得する引数の位置（0から始まる）。
:param defval: 引数が見つからない場合や変換できない場合に返すデフォルト値。
:returns: 変換された整数、またはデフォルト値。
"""
return pint(getarg(position, defval))
[ドキュメント]
def usage():
"""
プログラムの使用方法を標準出力に表示します。
:returns: None
"""
print("")
print("Usage: Variables in () are optional")
print(" python {}".format(sys.argv[0]))
print(" ex: python {} (basis)".format(sys.argv[0]))
print(" ex: python {} (band)".format(sys.argv[0]))
print(" ex: python {} (wf)".format(sys.argv[0]))
print("")
[ドキュメント]
def terminate(message = None):
"""
指定されたメッセージを表示し、プログラムを終了します。
終了前に`usage()`関数を呼び出して使用方法を表示します。
:param message: 表示するエラーメッセージ。
:returns: None
"""
if message is not None:
print("")
print(message)
print("")
usage()
print("")
exit()
[ドキュメント]
def make_list2(dim1, dim2):
"""
指定された次元の2次元リストを作成し、0.0で初期化します。
:param dim1: 1次元目のサイズ。
:param dim2: 2次元目のサイズ。
:returns: 0.0で初期化された2次元リスト。
"""
l = []
for i in range(dim1):
l.append([])
for j in range(dim2):
l[i].append(0.0)
return l
[ドキュメント]
def make_list3(dim1, dim2, dim3):
"""
指定された次元の3次元リストを作成し、0.0で初期化します。
:param dim1: 1次元目のサイズ。
:param dim2: 2次元目のサイズ。
:param dim3: 3次元目のサイズ。
:returns: 0.0で初期化された3次元リスト。
"""
l = make_list2(dim1, dim2)
for i1 in range(dim1):
for i2 in range(dim2):
l[i1][i2] = []
for i3 in range(dim3):
l[i1][i2].append(0.0)
return l
[ドキュメント]
def make_list4(dim1, dim2, dim3, dim4):
"""
指定された次元の4次元リストを作成し、0.0で初期化します。
:param dim1: 1次元目のサイズ。
:param dim2: 2次元目のサイズ。
:param dim3: 3次元目のサイズ。
:param dim4: 4次元目のサイズ。
:returns: 0.0で初期化された4次元リスト。
"""
l = make_list3(dim1, dim2, dim3)
for i1 in range(dim1):
for i2 in range(dim2):
for i3 in range(dim3):
l[i1][i2][i3] = []
for i4 in range(dim4):
l[i1][i2][i3].append(0.0)
return l
[ドキュメント]
def reduce(x, x0):
"""
xをx0で割った剰余を計算し、指定された範囲に減らします。
xが負の場合でも適切に動作します（Pythonの `%` 演算子とは異なる挙動）。
:param x: 減らす対象の数値。
:param x0: 周期の幅。
:returns: x0の範囲に減らされたxの値。
"""
n = int(x / x0)
if x < 0.0:
n -= 1
return x - x0 * n
#==============================================
# update default values by startup arguments
#==============================================
argv = sys.argv
#if len(argv) == 1:
# terminate()
mode = getarg (1, mode)
[ドキュメント]
def atominf(i, atoms):
"""
指定されたインデックスの原子情報を取り出します。
:param i: 原子リスト内でのインデックス。
:param atoms: 原子情報を格納したリスト。
:returns: 原子の名前、井戸の幅、基底ポテンシャル、有効質量。
"""
atom = atoms[i]
name = atom["name"]
wwidth = atom["wellwidth"]
Vbase = atom["Vbase"]
meff = atom["meff"]
return name, wwidth, Vbase, meff
[ドキュメント]
def siteinf(i, crystal):
"""
指定されたインデックスのサイト（結晶内の原子位置）情報を詳細に取り出します。
:param i: サイトリスト内でのインデックス。
:param crystal: サイト情報を格納したリスト。
:returns: サイト名、原子インデックス、原子辞書、サイトのX座標、井戸の幅、井戸の最小X座標、井戸の最大X座標、基底ポテンシャル、有効質量。
"""
site = crystal[i]
iatom = site["iatom"]
atom = atoms[iatom]
name = site["name"]
x = site["x0"]
wwidth = atom["wellwidth"]
Vbase = atom["Vbase"]
meff = atom["meff"]
return name, iatom, atom, x, wwidth, x - wwidth / 2.0, x + wwidth / 2.0, Vbase, meff
[ドキュメント]
def pot(x, Vvac, crystal):
"""
指定された位置におけるポテンシャルエネルギーを計算します。
x座標を単位格子aで周期的に還元し、結晶構造に基づいてポテンシャル値を決定します。
:param x: ポテンシャルを評価するX座標。
:param Vvac: 真空準位のポテンシャルエネルギー。
:param crystal: 結晶サイトの定義リスト。
:returns: X座標におけるポテンシャルエネルギー。
"""
xred = reduce(x, a)
nsites = len(crystal)
for i in range(nsites):
name, iatom, atom, x0, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
if x0 - wwidth/2.0 <= xred < x0 + wwidth/2.0:
return Vbase
else:
return Vvac
[ドキュメント]
def pot_atom(dx, Vvac, crystal):
"""
局所座標dxにおける原子ごとのポテンシャルエネルギーを計算します。
主にタイトバインディング計算の積分範囲で、原子を中心とした局所ポテンシャルを評価するために使用されます。
:param dx: 原子中心からの相対X座標。
:param Vvac: 真空準位のポテンシャルエネルギー。
:param crystal: 結晶サイトの定義リスト。
:returns: 相対X座標におけるポテンシャルエネルギー。
"""
nsites = len(crystal)
for i in range(nsites):
name, iatom, atom, x0, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
a = wwidth / 2.0
if -a <= dx < a:
return Vbase
else:
return Vvac
[ドキュメント]
def build_potential(xmin, xstep, n, Vvac, crystal):
"""
指定された範囲とステップでポテンシャルエネルギーのプロファイルを作成します。
:param xmin: X座標の最小値。
:param xstep: X座標のステップ幅。
:param n: サンプル点数。
:param Vvac: 真空準位のポテンシャルエネルギー。
:param crystal: 結晶サイトの定義リスト。
:returns: X座標の配列と対応するポテンシャルエネルギーの配列。
"""
xpot = np.empty(n)
ypot = np.empty(n)
for i in range(n):
xx = xmin + i * xstep
xpot[i] = xx
ypot[i] = pot(xx, Vvac, crystal)
return xpot, ypot
[ドキュメント]
def quantization_level_inifinite(Vbase, meff, wwidth, iL):
"""
無限に深い井戸型ポテンシャルにおける量子化されたエネルギー準位を計算します。
:param Vbase: 井戸内の基底ポテンシャルエネルギー。
:param meff: 有効質量（電子質量の倍数）。
:param wwidth: 井戸の幅。
:param iL: 量子数（1から始まる）。
:returns: 計算されたエネルギー準位。
"""
k = pi / wwidth
return Vbase + KE / meff * (k * iL)**2
[ドキュメント]
def cal_delta_even(E, w, V0, meff):
"""
有限量子井戸の偶関数解の特性方程式の左辺（デルタ関数）を計算します。
ポテンシャル井戸の幅`w`、井戸の深さ`V0`、有効質量`meff`が与えられたときの、エネルギー`E`に対応する特性方程式の一部を返します。
:param E: 井戸底からのエネルギー。
:param w: 井戸の幅。
:param V0: 井戸の深さ（真空準位と井戸底の差）。
:param meff: 有効質量。
:returns: デルタ関数（特性方程式の一部）の値。
"""
kb = sqrt(2.0 * meff * me * (V0 - E) * e) / hbar * 1.0e-10 # A^-1
kw = sqrt(2.0 * meff * me * E * e) / hbar * 1.0e-10 # A^-1
# kb = 0.0
delta = kb - kw * tan(kw * w / 2.0)
return delta
[ドキュメント]
def cal_delta_odd(E, w, V0, meff):
"""
有限量子井戸の奇関数解の特性方程式の左辺（デルタ関数）を計算します。
ポテンシャル井戸の幅`w`、井戸の深さ`V0`、有効質量`meff`が与えられたときの、エネルギー`E`に対応する特性方程式の一部を返します。
:param E: 井戸底からのエネルギー。
:param w: 井戸の幅。
:param V0: 井戸の深さ（真空準位と井戸底の差）。
:param meff: 有効質量。
:returns: デルタ関数（特性方程式の一部）の値。
"""
kb = sqrt(2.0 * meff * me * (V0 - E) * e) / hbar * 1.0e-10 # A^-1
kw = sqrt(2.0 * meff * me * E * e) / hbar * 1.0e-10 # A^-1
# kb = 0.0
delta = kb + kw / tan(kw * w / 2.0)
return delta
[ドキュメント]
def refine_E(cal_delta, E0, E1, nmaxiter, eps, w, V0, meff, IsPrint = 0):
"""
セカント法（割線法）を用いて、特性方程式の根（エネルギー準位）を精密化します。
`cal_delta`関数が符号を反転する区間`[E0, E1]`内で根を探索し、指定された精度`eps`または最大反復回数`nmaxiter`に達するまで反復します。
:param cal_delta: 根を探索するデルタ関数（偶関数または奇関数）。
:param E0: 探索範囲の開始エネルギー。
:param E1: 探索範囲の終了エネルギー。
:param nmaxiter: 最大反復回数。
:param eps: 収束判定の許容誤差。
:param w: 井戸の幅。
:param V0: 井戸の深さ。
:param meff: 有効質量。
:param IsPrint: 途中経過を表示するかどうかのフラグ（0で非表示）。
:returns: 精密化されたエネルギー準位、収束時のエネルギー変化、収束時のデルタ関数の値、または収束しなかった場合はNone。
"""
delta0 = cal_delta(E0, w, V0, meff)
# print(" 0:", E0, delta0)
delta1 = cal_delta(E1, w, V0, meff)
# print(" 1:", E1, delta1)
for i in range(nmaxiter):
if delta0 * delta1 > 0.0:
print("Error: delta0 ({:8.4f}, {:10.4g}) and delta1 ({:8.4f}, {:10.4g}) must have opposite signs"
.format(E0, delta0, E1, delta1))
exit()
E2 = (E0 + E1) / 2.0
delta2 = cal_delta(E2, w, V0, meff)
# print(" {}: E={:8.4f} - {:8.4f} => {:8.4f}".format(i, E0, E1, E2))
# print(" delta: {:12.4g} - {:12.4g} => {:12.4g}".format(delta0, delta1, delta2))
if abs(E2 - E0) < eps:
if IsPrint:
print(" converged at E = {:12.6g} with dE = {:12.6g} delta = {:12.6g}"
.format(E2, E2 - E0, delta2))
return E2, E2 - E0, delta2
else:
if delta0 * delta2 > 0.0:
E0 = E2
delta0 = delta2
else:
E1 = E2
delta1 = delta2
continue
else:
print(" Not converged for {} iterations.".format(nmaxiter))
print(" E = {:12.6g} with dE = {:12.6g} delta = {:12.6g}".format(E2, dE, delta2))
return None, None, None
[ドキュメント]
def find_eigenstates_base(parity, Vvac, Vbase, wwidth, meff, Emin, Emax, nE, nmaxiter, eps, IsPrint = 0):
"""
指定されたパリティ（偶奇性）を持つ有限量子井戸の固有状態エネルギーを探索します。
`Emin`から`Emax`までのエネルギー範囲を`nE`点で走査し、`cal_delta`関数の符号反転を検出して、`refine_E`を用いて固有エネルギーを精密化します。
:param parity: 固有状態のパリティ (0:偶関数, 1:奇関数)。
:param Vvac: 真空準位のポテンシャルエネルギー。
:param Vbase: 井戸内の基底ポテンシャルエネルギー。
:param wwidth: 井戸の幅。
:param meff: 有効質量。
:param Emin: エネルギー探索範囲の最小値。
:param Emax: エネルギー探索範囲の最大値。
:param nE: エネルギー探索の分割数。
:param nmaxiter: セカント法の最大反復回数。
:param eps: セカント法の収束判定の許容誤差。
:param IsPrint: 途中経過を表示するかどうかのフラグ（0で非表示）。
:returns: 見つかった固有状態のリスト（各要素はエネルギーとパリティを含む辞書）。
"""
V0 = Vvac - Vbase
Estep = (Emax - Emin) / (nE - 1)
if parity:
cal_delta = cal_delta_odd
else:
cal_delta = cal_delta_even
d0 = None
iband = 0
orbitals = []
for iE in range(nE):
E = Emin + iE * Estep
if E <= Vbase:
continue
if Vvac <= E:
break
delta = cal_delta(E - Vbase, wwidth, V0, meff)
if d0 is None:
d0 = delta
continue
# cal_delta function flips its sign when delta => 0 and delta => infinite
# the latter case must be removed
if d0 * delta < 0.0 and delta < 10.0:
# print("0: ", E-Estep, d0)
# print("1: ", E, delta)
E0, dE, delta0 = refine_E(cal_delta, E - Vbase - Estep, E - Vbase, nmaxiter, eps, wwidth, V0, meff, IsPrint = 0)
E = E0 + Vbase
if IsPrint:
print(" E[{}]={:12.6g} eV dE={:12.6g} delta={:12.6g}".format(iband, E, dE, delta0))
orbitals.append({"parity": parity, "E": E})
d0 = delta
E += Estep
iband += 1
if d0 * delta < 0.0:
d0 = delta
return orbitals
[ドキュメント]
def find_eigenstates(Vvac, Vbase, wwidth, meff, Emin, Emax, nE, nmaxiter, eps):
"""
有限量子井戸の偶関数および奇関数の固有状態エネルギーとその波動関数パラメータを計算します。
`find_eigenstates_base`を使用して偶関数と奇関数の固有エネルギーを探索し、それらをエネルギー順にソートします。
各固有状態に対して、波動関数の係数(`Cb`, `Cw`)と波長(`kb`, `kw`)を計算し、結果に含めます。
:param Vvac: 真空準位のポテンシャルエネルギー。
:param Vbase: 井戸内の基底ポテンシャルエネルギー。
:param wwidth: 井戸の幅。
:param meff: 有効質量。
:param Emin: エネルギー探索範囲の最小値。
:param Emax: エネルギー探索範囲の最大値。
:param nE: エネルギー探索の分割数。
:param nmaxiter: セカント法の最大反復回数。
:param eps: セカント法の収束判定の許容誤差。
:returns: 固有状態のリスト（各要素はエネルギー、パリティ、波動関数パラメータを含む辞書）。
"""
a = wwidth / 2.0
orbitals_o = find_eigenstates_base(1, Vvac, Vbase, wwidth, meff, Emin, Emax, nE, nmaxiter, eps)
orbitals_e = find_eigenstates_base(0, Vvac, Vbase, wwidth, meff, Emin, Emax, nE, nmaxiter, eps)
orbitals = sorted(orbitals_o + orbitals_e, key = lambda x: x["E"])
norbitals = len(orbitals)
for i in range(norbitals):
p = orbitals[i]["parity"]
E = orbitals[i]["E"]
# print(" {:3d}: parity={:4} E={:12.6g} eV".format(i, parity_str[p], E))
V0 = Vvac - Vbase
kb = sqrt(2.0 * meff * me * (V0 - E) * e) / hbar * 1.0e-10 # A^-1
kw = sqrt(2.0 * meff * me * E * e) / hbar * 1.0e-10 # A^-1
Cw = 1.0 / sqrt(wwidth/2.0 + 1.0 / kb)
if p == 0:
Cb = Cw * cos(kw * a)
else:
Cb = Cw * sin(kw * a)
orbitals[i]["kb"] = kb
orbitals[i]["kw"] = kw
orbitals[i]["Cb"] = Cb
orbitals[i]["Cw"] = Cw
orbitals[i]["Vvac"] = Vvac
orbitals[i]["Vbase"] = Vbase
orbitals[i]["wwidth"] = wwidth
return orbitals
[ドキュメント]
def basis(dx, orbital, Nx = 0):
"""
指定された軌道（固有状態）の波動関数値を計算します。
有限量子井戸ポテンシャルにおける波動関数の解析解に基づき、中心からの相対距離`dx`における波動関数の値を返します。
:param dx: 井戸の中心からの相対X座標。
:param orbital: 固有状態のパラメータ（パリティ、幅、エネルギー、係数など）を含む辞書。
:param Nx: この引数は現在の実装では使用されていません。
:returns: 波動関数値。
"""
p = orbital["parity"]
wwidth = orbital["wwidth"]
a = wwidth / 2.0
E = orbital["E"]
kb = orbital["kb"]
kw = orbital["kw"]
Cb = orbital["Cb"]
Cw = orbital["Cw"]
if p == 0:
if dx < -a:
return Cb * exp(kb * (dx + a))
elif dx > a:
return Cb * exp(-kb * (dx - a))
else:
return Cw * cos(kw * dx)
else:
if dx < -a:
return -Cb * exp(kb * (dx + a))
elif dx > a:
return Cb * exp(-kb * (dx - a))
else:
return Cw * sin(kw * dx)
[ドキュメント]
def cal_basis(x0, wfmin, wfmax, nwf, orbital):
"""
特定の軌道（固有状態）の波動関数を、指定されたX範囲で計算します。
`wfmin`から`wfmax`までの範囲を`nwf`点でサンプリングし、各点で`basis`関数を呼び出して波動関数プロファイルを作成します。
結果のX座標は`x0`でオフセットされます。
:param x0: 波動関数の中心X座標。
:param wfmin: 計算範囲の相対X座標最小値。
:param wfmax: 計算範囲の相対X座標最大値。
:param nwf: 計算点数。
:param orbital: 固有状態のパラメータを含む辞書。
:returns: 絶対X座標のリストと対応する波動関数値のリスト。
"""
wfstep = (wfmax - wfmin) / (nwf - 1)
xx = []
ywf = []
# print(" E={} eV parity={}".format(Edic["E"], Edic["parity"]))
# print(" x={} - {}".format(wfmin, wfmax))
for i in range(nwf):
x = wfmin + i * wfstep
wf = basis(x, orbital)
xx.append(x0 + x)
ywf.append(wf)
# print(x, wf)
return xx, ywf
[ドキュメント]
def cal_cdiff2(h, f):
"""
複素数値配列の2階中心差分を計算します。
空間間隔`h`と複素数値関数`f`の配列が与えられたとき、各点での2階微分（ラプラシアン）を近似します。
境界では隣接する点を使用します。
:param h: 空間ステップ幅。
:param f: 複素数値関数の配列。
:returns: 2階微分値の複素数値配列。
"""
n = len(f)
diff2 = np.empty(n, dtype = complex)
for i in range(n):
if i == 0:
idx = 1
elif i == n - 1:
idx = n - 2
else:
idx = i
d2r = (f[idx+1].real - 2.0 * f[idx].real + f[idx-1].real) / h / h
d2i = (f[idx+1].imag - 2.0 * f[idx].imag + f[idx-1].imag) / h / h
diff2[i] = d2r + 1.0j * d2i
return diff2
[ドキュメント]
def integrate2(h, y1, y2):
"""
2つの関数`y1`と`y2`の積の数値積分を台形近似で計算します。
`h`をステップ幅として、`y1[i] * y2[i]`の値を合計し、`h`を乗じて積分値を近似します。
:param h: 積分ステップ幅。
:param y1: 1つ目の関数の値のリスト。
:param y2: 2つ目の関数の値のリスト。
:returns: 積分値。
"""
S = 0.0
for i in range(len(y1)-1):
S += y1[i] * y2[i]
return S * h
[ドキュメント]
def integrate3(h, y1, y2, y3):
"""
3つの関数`y1`, `y2`, `y3`の積の数値積分を台形近似で計算します。
`h`をステップ幅として、`y1[i] * y2[i] * y3[i]`の値を合計し、`h`を乗じて積分値を近似します。
:param h: 積分ステップ幅。
:param y1: 1つ目の関数の値のリスト。
:param y2: 2つ目の関数の値のリスト。
:param y3: 3つ目の関数の値のリスト。
:returns: 積分値。
"""
S = 0.0
for i in range(len(y1)-1):
S += y1[i] * y2[i] * y3[i]
return S * h
[ドキュメント]
def cal_E(orbital, xmin, xmax, nx):
"""
波動関数とポテンシャルから、固有状態の運動エネルギーとポテンシャルエネルギーを数値的に計算し、全エネルギーを導出します。
与えられた軌道（固有状態）の波動関数を使用し、指定されたX範囲でポテンシャルエネルギーと波動関数の2階微分を計算します。
これらを用いて、規格化積分、運動エネルギー、ポテンシャルエネルギー、および全エネルギーを数値積分によって計算します。
:param orbital: 計算対象の固有状態のパラメータを含む辞書。
:param xmin: 積分範囲の最小X座標。
:param xmax: 積分範囲の最大X座標。
:param nx: 積分点数。
:returns: 規格化積分値、運動エネルギー、ポテンシャルエネルギー、全エネルギー。
"""
global Vvac, crystal
xstep = (xmax - xmin) / (nx - 1)
xx = [xmin + i * xstep for i in range(nx)]
ypot = np.empty(nx, dtype = complex)
ywf = np.empty(nx, dtype = complex)
for i in range(nx):
x = xx[i]
ypot[i] = pot_atom(x, Vvac, crystal)
ywf[i] = basis(x, orbital)
ydiff2 = cal_cdiff2(xstep, ywf)
"""
print("")
print("x, U(x), wf(x), d2wf/dx2")
for i in range(nx):
print("{:8.4f}\t{:8.4f}\t{:8.4f}+j{:8.4f}\t{:8.4f}+j{:8.4f}"
.format(xx[i], ypot[i], ywf[i].real, ywf[i].imag, ydiff2[i].real, ydiff2[i].imag))
print("")
"""
Sii = integrate2(xstep, ywf.conjugate(), ywf).real
K = -KE * integrate2(xstep, ywf.conjugate(), ydiff2).real
U = integrate3(xstep, ywf.conjugate(), ypot, ywf).real
E = K + U
# print("K=", Sii, K, U, E)
# exit()
return Sii, K, U, E
[ドキュメント]
def cal_E_analytical(orbital):
"""
有限量子井戸の固有状態に対して、規格化積分、運動エネルギー、ポテンシャルエネルギー、全エネルギーを解析的に計算します。
井戸内部（well）と障壁外部（barrier）に分けて解析的に積分を実行し、それぞれの寄与を計算します。
:param orbital: 計算対象の固有状態のパラメータを含む辞書。
:returns: 障壁部のSii、井戸部のSii、全Sii、障壁部のK、井戸部のK、全K、障壁部のU、井戸部のU、全U、障壁部のE、井戸部のE、全E。
"""
p = orbital["parity"]
wwidth = orbital["wwidth"]
a = wwidth / 2.0
E = orbital["E"]
kb = orbital["kb"]
kw = orbital["kw"]
Cb = orbital["Cb"]
Cw = orbital["Cw"]
Vvac = orbital["Vvac"]
Vbase = orbital["Vbase"]
# x > a: Cb * exp(kb*a) * exp(-kb * x)
Siib = Cb * Cb / kb * exp(2.0 * kb * a) * exp(-2.0 * kb * a)
if p == 0:
# -a < x < a: Cw * cos(kw*x)
# integ = integ(C*C * cos^2(kx)) = C*C*0.5*(cos(2kx)+1)]=C*C*0.5*(0.5/k*sin(2kx)+x)_[-a,a]
# = 0.5 * C*C * [0.5*(sin(2ka)-sin(-2ka)) / k + 2a]
Siiw = Cw * Cw * 0.5 * (0.5*(sin(2.0 * kw * a)-sin(-2.0 * kw * a)) / kw + 2.0 * a)
else:
# -a < x < a: Cw * sin(kw * dx)
# integ = integ(C*C * sin^2(kx)) = C*C*0.5*(1-cos(2kx))]=C*C*0.5*(-0.5/k*sin(2kx)+x)_[-a,a]
# = 0.5 * C*C * [-0.5*(sin(2.0 * kw * a)-sin(-2.0 * kw * a)) / kw + 2.0 * a]
Siiw = Cw * Cw * 0.5 * (-0.5*(sin(2.0 * kw * a)-sin(-2.0 * kw * a)) / kw + 2.0 * a)
# x > a: d2f/dx2 = Cb * exp(kb*a) * kb*kb * exp(-kb * x) = kb*kb * f
Kiib = -KE * kb * kb * Siib
# -a < x < a: d2f/dx2 = -Cw * kw*kw * cos(kw*x) = -kw*kw * f
Kiiw = KE * kw * kw * Siiw
Ub = Vvac * Siib
Uw = Vbase * Siiw
return Siib, Siiw, Siib+Siiw, Kiib, Kiiw, Kiib+Kiiw, Ub, Uw, Ub+Uw, Kiib+Ub, Kiiw+Uw, Kiib+Ub+Kiiw+Uw
[ドキュメント]
def wf_innerproduct(x, y1, y2):
"""
2つの波動関数`y1`と`y2`の内積を数値的に計算します。
与えられたX座標`x`と関数値`y1`, `y2`を用いて、台形近似により内積を計算します。
:param x: X座標のリスト。
:param y1: 1つ目の波動関数の値のリスト。
:param y2: 2つ目の波動関数の値のリスト。
:returns: 2つの波動関数の内積。
"""
S = 0.0
for i in range(len(x)-1):
S += (y1[i].conjugate() * y2[i] + y1[i+1].conjugate() * y2[i+1]) * (x[i+1] - x[i])
S *= 0.5
return S
[ドキュメント]
def cal_tij(xmin, xmax, nx, orbital1, orbital2, x01, x02, Nx):
"""
2つの異なるサイトにある波動関数間の遷移積分（ホッピング積分）t_ijを計算します。
サイト1の波動関数1（中心x01）とサイト2の波動関数2（中心x02 + Nx*a）がオーバーラップする領域で、
ポテンシャルを介した積を数値積分します。
:param xmin: 積分範囲の最小X座標。
:param xmax: 積分範囲の最大X座標。
:param nx: 積分点数。
:param orbital1: 1つ目のサイトの固有状態のパラメータを含む辞書。
:param orbital2: 2つ目のサイトの固有状態のパラメータを含む辞書。
:param x01: 1つ目のサイトのX座標。
:param x02: 2つ目のサイトのX座標。
:param Nx: サイト2がサイト1からどれだけ周期的に離れているかを示す整数（Nx*a）。
:returns: 計算された遷移積分値（複素数）。
"""
global Vvac, crystal, a
if orbital1 is orbital2 and Nx == 0:
return orbital1["E"]
# print("orb {}:{} - {}:{}".format(orbital1["E"], orbital1["parity"], orbital2["E"], orbital2["parity"]))
xstep = (xmax - xmin) / (nx - 1)
S = 0.0
for i in range(nx):
x = xmin + i * xstep
V = pot(x, Vvac, crystal)
b1 = basis(x - x01, orbital1).conjugate()
b2 = basis(x - x02 - Nx*a, orbital2)
S += b1 * V * b2
# if abs(b1) > 1.0e-3 and abs(b2) > 1.0e-3 and V > 0.0:
# print(" {}({:8.4f}):{:8.4f}, {:8.4f}, {:8.4f} S={:8.4f}".format(i, x, b1, V, b2, S))
S *= xstep
return S
[ドキュメント]
def cal_Sij(xmin, xmax, nx, orbital1, orbital2, x01, x02, Nx):
"""
2つの異なるサイトにある波動関数間のオーバーラップ積分S_ijを計算します。
サイト1の波動関数1（中心x01）とサイト2の波動関数2（中心x02 + Nx*a）がオーバーラップする領域で、
積を数値積分します。
:param xmin: 積分範囲の最小X座標。
:param xmax: 積分範囲の最大X座標。
:param nx: 積分点数。
:param orbital1: 1つ目のサイトの固有状態のパラメータを含む辞書。
:param orbital2: 2つ目のサイトの固有状態のパラメータを含む辞書。
:param x01: 1つ目のサイトのX座標。
:param x02: 2つ目のサイトのX座標。
:param Nx: サイト2がサイト1からどれだけ周期的に離れているかを示す整数（Nx*a）。
:returns: 計算されたオーバーラップ積分値（複素数）。
"""
global crystal, a
if orbital1 is orbital2 and Nx == 0:
return 1.0
xstep = (xmax - xmin) / (nx - 1)
S = 0.0
for i in range(nx):
x = xmin + i * xstep
S += basis(x - x01, orbital1).conjugate() * basis(x - x02 - Nx*a, orbital2)
S *= xstep
return S
[ドキュメント]
def build_basisset(crystal):
"""
結晶構造に基づいて、タイトバインディング法の基底関数セットを構築します。
各サイトとそこに含まれる原子の各軌道（固有状態）を組み合わせ、一意の基底状態としてリストにまとめます。
:param crystal: 結晶サイトの定義リスト。
:returns: 基底関数のリスト（各要素はサイト、軌道、原子情報、座標、固有状態パラメータを含む辞書）。
"""
nsites = len(crystal)
b = []
for isite in range(nsites):
name, iatom, atom, x0, wwidth, xwmin, xwmax, Vbase, meff = siteinf(isite, crystal)
orbitals = atom["orbitals"]
norbitals = len(orbitals)
for io in range(norbitals):
b.append({"isite": isite, "iorb": io, "atom": atom, "name": name, "x0": x0, "orbital": orbitals[io]})
return b
[ドキュメント]
def search_neighbors(rmax, Nxmax, crystal):
"""
指定された範囲内で、各サイトの近接サイトを検索します。
各サイトについて、格子周期`a`を用いて`Nxmax`周期分までの範囲で他のサイトを探索し、`rmax`以内に存在するサイトを近接サイトとしてリストアップします。
:param rmax: 近接サイトと見なす最大距離。
:param Nxmax: 探索する格子周期の最大範囲（-Nxmax から +Nxmax）。
:param crystal: 結晶サイトの定義リスト。
:returns: 各サイトの近接サイトリストのリスト。
"""
global a
nsites = len(crystal)
nlist = []
for isite1 in range(nsites):
nlist.append([])
name1, iatom1, atom1, x01, wwidth, xwmin, xwmax, Vbase, meff = siteinf(isite1, crystal)
for isite2 in range(nsites):
name2, iatom2, atom2, x02, wwidth, xwmin, xwmax, Vbase, meff = siteinf(isite2, crystal)
for inx in range(-Nxmax, +Nxmax+1):
r = abs(x02 + a * inx - x01)
# if 1.0e-6 < r and r < rmax:
if r < rmax:
nlist[isite1].append({"isite2": isite2, "Nx": inx,
"atom2": atom2, "name2": name2, "x02": x02, "r": r})
return nlist
[ドキュメント]
def solve_eigenstates(k, basisset, neighborlist, IsPrint = 0):
"""
指定された波数`k`に対して、タイトバインディング法によるハミルトニアンとオーバーラップ行列を構築し、固有状態（バンドエネルギー）を解きます。
構築された基底セットと近接サイト情報を用いて、ハミルトニアン行列`Fij`とオーバーラップ行列`Sij`を計算します。
その後、一般化固有値問題 `Fij * c = E * Sij * c` を解く代わりに、`Sij`の逆行列を`Fij`に左から掛けた行列の固有値を計算します。
:param k: 波数（π/a単位）。
:param basisset: 基底関数のリスト。
:param neighborlist: 各サイトの近接サイトのリスト。
:param IsPrint: 詳細な計算過程を表示するかどうかのフラグ（0で非表示）。
:returns: 計算されたバンドエネルギーの配列（実部）、固有ベクトルの配列（現在の実装では空）、および計算結果の詳細を含む辞書。
"""
global tijxmin, tijxmax, ntijx
tijxstep = (tijxmax - tijxmin) / (ntijx - 1)
# hij = make_list4(nsites, nbmax, nsites, nbmax)
nbs = len(basisset)
Fij = np.zeros([nbs, nbs], dtype = complex)
Sij = np.zeros([nbs, nbs], dtype = complex)
# diagonal part
# for ib1 in range(nbs):
# bs1 = basisset[ib1]
# orbital1 = bs1["orbital"]
# E1 = orbital1["E"]
# Fij[ib1][ib1] = E1
# Sij[ib1][ib1] = 1.0
# non-diagonal part
for ib1 in range(nbs):
bs1 = basisset[ib1]
isite1 = bs1["isite"]
iorb1 = bs1["iorb"]
name1 = bs1["name"]
atom1 = bs1["atom"]
x01 = bs1["x0"]
orbital1 = bs1["orbital"]
E1 = orbital1["E"]
nbl1 = neighborlist[isite1]
for ib2 in range(ib1, nbs):
# for ib2 in range(ib1+1, nbs):
bs2 = basisset[ib2]
isite2 = bs2["isite"]
iorb2 = bs2["iorb"]
name2 = bs2["name"]
atom2 = bs2["atom"]
x02 = bs2["x0"]
orbital2 = bs2["orbital"]
E2 = orbital2["E"]
# if ib1 == ib2:
# Fij[ib1][ib2] = E1
# Sij[ib1][ib2] = 1.0
for inb in range(len(nbl1)):
nb = nbl1[inb]
if isite2 != nb["isite2"]:
continue
# if ib1 != ib2 and nb["r"] > 0.1:
# pass
# else:
# continue
Nx = nb["Nx"]
if IsPrint >= 1:
print(" base {}-{} {}[L{}]-{}[L{}]({:+d}a)): "
.format(ib1, ib2, name1, iorb1, name2, iorb2, Nx), end = '')
print( " ({:6.4g}eV)".format(E1), end = '')
print(" - ({:6.4g}eV)".format(E2))
# ximin = max(x01 + tijxmin, x02 + tijxmin)
ximin = min(x01 + tijxmin, x02 + tijxmin)
# ximax = min(x01 + tijxmax, x02 + tijxmax)
ximax = max(x01 + tijxmax, x02 + tijxmax)
nxi = int((ximax - ximin) / tijxstep + 1.00001)
if IsPrint >= 2:
print(" tij integration range: {:6.4g} - {:6.4g}, {:6.4g} step, {} points"
.format(ximin, ximax, tijxstep, nxi))
tij_ = cal_tij(ximin, ximax, nxi, orbital1, orbital2, x01, x02, Nx)
sij_ = cal_Sij(ximin, ximax, nxi, orbital1, orbital2, x01, x02, Nx)
phase = exp(1.0j * k * pi2/a * (x02 + Nx*a - x01))
Fij[ib1][ib2] += tij_ * phase
Sij[ib1][ib2] += sij_ * phase
if IsPrint >= 2:
print(" phase=", phase, Nx)
print(" tij = {:12.6g}".format(tij_))
print(" sij = {:12.6g}".format(sij_))
for ib1 in range(nbs):
for ib2 in range(ib1+1, nbs):
Fij[ib2][ib1] = Fij[ib1][ib2].conjugate()
Sij[ib2][ib1] = Sij[ib1][ib2].conjugate()
if IsPrint >= 1:
print("")
print("Fij=\n", Fij)
print("Sij=\n", Sij)
Fij = LA.inv(Sij) @ Fij
ei = LA.eigvals(Fij)
ci = []
# ei, ci = LA.eigh(Fij)
if IsPrint >= 2:
print("w=", w)
print("v=", v)
res = {"ei": ei, "ci": ci, "k": k, "Fij": Fij, "Sij": Sij}
return ei.real, ci, res
[ドキュメント]
def wf():
"""
波動関数プロットモードの処理を実行します。
現在、この関数は実装されておらず、`terminate()`を呼び出して終了します。
:returns: None
"""
global mode
global a
global V0, crystal
global xmin, xmax, nx
global tijxmin, tijxmax, ntijx
terminate()
[ドキュメント]
def band():
"""
バンド構造計算モードの処理を実行します。
結晶内の有限量子井戸の固有状態を計算し、それらを基底関数として使用して、
タイトバインディング法によるバンド構造を計算します。
計算されたバンド構造とポテンシャル、井戸内準位をプロットします。
:returns: None
"""
global mode
global a
global V0, crystal
global xmin, xmax, nx
global tijxmin, tijxmax, ntijx
nsites = len(crystal)
tijxstep = (tijxmax - tijxmin) / (ntijx - 1)
xmin = min(tijxmin, a + tijxmin)
xmax = max(tijxmax, a + tijxmax)
xstep = (xmax - xmin) / (nx - 1)
print("")
print("=== Input parameterss ===")
print("mode:", mode)
print("Crystal crystal")
print(" a=", a, "A")
print(" Vacuum potential: {:10.6f} eV".format(Vvac))
for i in range(nsites):
name, iatom, atom, x, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
print(" {:8s}: x={:10.6f} A".format(name, x))
print(" well: V(x)={:10.6f} eV from {:10.6f} to {:10.6f} A"
.format(Vbase, xwmin, xwmax))
print("X range to calculate tij")
print(" {:10.6f} - {:10.6f} A at {:10.6f} A step, {} points".format(tijxmin, tijxmax, tijxstep, ntijx))
xpot, ypot = build_potential(xmin, xstep, nx, Vvac, crystal)
print("")
print("=== for finite potential ===")
print(" Vacuum potential: {:10.6f} eV".format(Vvac))
atom_orbitals = []
xWFlist = []
yWFlist = []
for i in range(nsites):
name, iatom, atom, x0, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
print(" Atom {}: {}".format(i, name))
orbitals = find_eigenstates(Vvac, Vbase, wwidth, meff, Emin, Emax, nE, nmaxiter, eps)
WFlist_ = []
for ie in range(len(orbitals)):
orbital = orbitals[ie]
p = orbital["parity"]
E = orbital["E"]
print(" Level {:3d}: parity={:4} E={:12.6g} eV".format(i, parity_str[p], E))
xwf, ywf = cal_basis(x0, tijxmin, tijxmax, ntijx, orbital)
wf2 = wf_innerproduct(xwf, ywf, ywf)
orbital["wf2"] = wf2
print(" wf2={}".format(wf2))
WFlist_.append(ywf)
atom["orbitals"] = orbitals
atom["xwf"] = xwf
atom["wfs"] = WFlist_
# atom_orbitals.append(orbitals)
xWFlist.append(xwf)
yWFlist.append(WFlist_)
print("")
print("=== Neighbor sites to calculate tij ===")
neighborlist = search_neighbors(rmax_neighbor, Nxmax_neighbor, crystal)
nbmax = max([len(nb) for nb in neighborlist])
print("nbmax=", nbmax)
print("rmax: {:8.4g} A".format(rmax_neighbor))
print("search Nx range: {} - {}".format(-Nxmax_neighbor, Nxmax_neighbor))
for isite1 in range(nsites):
name1, iatom1, atom1, x01, wwidth, xwmin, xwmax, Vbase, meff = siteinf(isite1, crystal)
print(" site {:03d} {:4}: x=({:8.4g})".format(isite1, name1, x01))
nbl = neighborlist[isite1]
for inb in range(len(nbl)):
nb = nbl[inb]
print(" - site {:03d} {:4}: x=({:8.4g} + {:+1d}*a) A r={:8.4g} A"
.format(nb["isite2"], nb["name2"], nb["x02"], nb["Nx"], nb["r"]))
print("")
print("=== Basis set ===")
basisset = build_basisset(crystal)
for i in range(len(basisset)):
b = basisset[i]
print(" {:03d}: site {:03d}:orb {:03d} {:4}".format(i, b["isite"], b["iorb"], b["name"]))
kstep = (kmax - kmin) / (nk - 1)
xk = [kmin + i * kstep for i in range(nk)]
yE0 = []
yE1 = []
yE = []
if nsites == 1:
print("")
print("=== Fij and non-hybridized bands for one-atom crystal ===")
print("k range: {} - {}, {} step (in pi/a), {} points".format(kmin, kmax, kstep, nk))
name1, iatom1, atom1, x01, wwidth, xwmin, xwmax, Vbase, meff = siteinf(0, crystal)
# ximin = max(x01 + tijxmin, x01 + tijxmin)
ximin = min(x01 + tijxmin, x01 + tijxmin)
# ximax = min(x01 + tijxmax, x01 + tijxmax)
ximax = max(x01 + tijxmax, x01 + tijxmax)
nxi = int((ximax - ximin) / tijxstep + 1.00001)
for i in range(nk):
k = xk[i]
phase = pi2 * k
phasem = exp(-1.0j * phase)
phasep = exp( 1.0j * phase)
ei0 = []
ei1 = []
for ib1 in range(len(basisset)):
bs1 = basisset[ib1]
iorb1 = bs1["iorb"]
print(" iorb=", iorb1)
orbital1 = bs1["orbital"]
tii = cal_tij(ximin, ximax, nxi, orbital1, orbital1, x01, x01, 0)
sii = cal_Sij(ximin, ximax, nxi, orbital1, orbital1, x01, x01, 0)
tijm = cal_tij(ximin, ximax, nxi, orbital1, orbital1, x01, x01, -1)
sijm = cal_Sij(ximin, ximax, nxi, orbital1, orbital1, x01, x01, -1)
tijp = cal_tij(ximin, ximax, nxi, orbital1, orbital1, x01, x01, +1)
sijp = cal_Sij(ximin, ximax, nxi, orbital1, orbital1, x01, x01, +1)
# Fij = phasem * tijm + tii + phasep * tijp
# Sij = phasem * sijm + sii + phasep * sijp
Fij = 2.0 * cos(phase) * tijp + tii
Sij = 2.0 * cos(phase) * sijp + sii
print(" at k={:8.4g}: Fij={:10.6g} Sij={:10.6g} Fij/Sij={:10.6g}"
.format(k, Fij, Sij, Fij / Sij))
# print(" phase=", phasem, phasep)
print(" diagonal: tii={:10.6g} sii={:10.6g}".format(tii, sii))
print(" to Nx=-1: tij={:10.6g} sij={:10.6g}".format(tijm, sijm))
print(" to Nx=+1: tij={:10.6g} sij={:10.6g}".format(tijp, sijp))
ei0.append(Fij)
ei1.append(Fij / Sij)
yE0.append(ei0)
yE1.append(ei1)
# exit()
print("")
print("=== Band ===")
print("k range: {} - {}, {} step (in pi/a), {} points".format(kmin, kmax, kstep, nk))
for i in range(nk):
k = xk[i]
ei, ci, res = solve_eigenstates(k, basisset, neighborlist, 0)
print("k: {:10.6f}".format(k), ":", ei)
yE.append(ei)
# print(" Fij=", res["Fij"].real)
# print(" Sij=", res["Sij"].real)
print("")
print("plot")
print("")
fig, ax = plt.subplots(1, 2)
ax1 = ax[0]
ax2 = ax[1]
# ax2をax1に関連させる
# ax2 = ax1.twinx()
ax1.plot(xpot, ypot, linewidth = 0.5) #, label = '$V$($x$)')
# Plot energy levels and wave functions for finite well potential
for i in range(nsites):
name, iatom, atom, x0, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
orbitals = atom["orbitals"]
# orbitals = atom_orbitals[i]
norbitals = len(orbitals)
for ie in range(norbitals):
orbital = orbitals[ie]
E = orbital["E"]
ax1.plot([xwmin, xwmax], [E, E], linestyle = 'dashed', color = 'b', linewidth = 1.0)
kwf = wfplotmax / max(yWFlist[i][ie])
ywf = [yWFlist[i][ie][j] * kwf + E for j in range(ntijx)]
ax1.plot(xWFlist[i], ywf, linestyle = '-', color = 'g', linewidth = 0.5)
ax2.set_xlim([-0.5, 0.5])
ax2.set_ylim(Erange)
ax2.plot(xk, yE, label = 'hybridized', linestyle = '', marker = 'o', markersize = 5.0,
markerfacecolor = 'none', markeredgecolor = 'black', markeredgewidth = 0.5)
if nsites == 1:
ax2.plot(xk, yE0, label = 'Fij', linestyle = '-', color = 'red', linewidth = 0.5)
ax2.plot(xk, yE1, label = 'Fij/Sij', linestyle = 'dashed', color = 'blue', linewidth = 0.5)
ax1.set_xlabel("x ($\AA$)", fontsize = fontsize)
ax1.set_ylabel("$V$($x$), $\phi$($x$)", fontsize = fontsize)
ax1.legend(fontsize = legend_fontsize)
ax1.tick_params(labelsize = fontsize)
ax2.set_xlabel("$k$ $(\pi$$/a)$", fontsize = fontsize)
ax2.set_ylabel("E (eV)", fontsize = fontsize)
# ax2.legend(fontsize = legend_fontsize)
ax2.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
terminate()
[ドキュメント]
def view_basis():
"""
基底関数表示モードの処理を実行します。
設定された結晶内の原子ごとのポテンシャル、無限井戸・有限井戸の量子化準位、
および有限井戸の波動関数を計算し、プロットします。
また、解析的および数値的なエネルギー計算結果も表示します。
:returns: None
"""
global mode
global a
global V0, atoms, crystal
global xmin, xmax, nx
natoms = len(atoms)
nsites = len(crystal)
xstep = (xmax - xmin) / (nx - 1)
print("")
print("=== Input parameterss ===")
print("mode:", mode)
print("Atoms:")
for i in range(natoms):
name, wwidth, Vbase, meff = atominf(i, atoms)
print(" {:8s}: well: width={:10.6f} A Vbase={:10.6f} eV me={:10.6f}me"
.format(name, wwidth, Vbase, meff))
print("Crystal structure:")
print(" a=", a, "A")
print(" Vacuum potential: {:10.6f} eV".format(Vvac))
for i in range(nsites):
name, iatom, atom, x, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
print(" {:8s}: x={:10.6f} A".format(name, x))
print(" well: V(x)={:10.6f} eV from {:10.6f} to {:10.6f} A"
.format(Vbase, xwmin, xwmax))
print("X range to plot V(x) and basis")
print(" {:10.6f} - {:10.6f} A at {:10.6f} A step, {} points".format(xmin, xmax, xstep, nx))
xpot, ypot = build_potential(xmin, xstep, nx, Vvac, crystal)
print("")
print("=== for infinite potential ===")
for i in range(natoms):
name, wwidth, Vbase, meff = atominf(i, atoms)
print(" {:8s}: well: width={:10.6f} A Vbase={:10.6f} eV me={:10.6f}me"
.format(name, wwidth, Vbase, meff))
for iL in range(1, 100):
E = quantization_level_inifinite(Vbase, meff, wwidth, iL)
if E > Vvac:
break
print(" Level {:3d}: E={:10.6f} eV E-Vbase={:10.6f} eV".format(iL, E, E - Vbase))
print("")
print("=== for finite potential ===")
print(" Vacuum potential: {:10.6f} eV".format(Vvac))
# atom_orbitals = []
xWFlist = []
yWFlist = []
for i in range(nsites):
name, iatom, atom, x0, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
print(" {:8s}: well: width={:10.6f} A Vbase={:10.6f} eV me={:10.6f}me"
.format(name, wwidth, Vbase, meff))
xstep = (xmax - xmin) / (nx - 1)
print(" integration x range: {} - {}, {} step, n={}".format(xmin, xmax, xstep, nx))
WFlist_ = []
orbitals = find_eigenstates(Vvac, Vbase, wwidth, meff, Emin, Emax, nE, nmaxiter, eps)
norbitals = len(orbitals)
for ie in range(norbitals):
orbital = orbitals[ie]
p = orbital["parity"]
E = orbital["E"]
print(" Level {:3d}: parity={:4} E={:12.6g} eV".format(ie, parity_str[p], E))
xwf, ywf = cal_basis(x0, -wwidth/2.0 - wmargin, wwidth/2.0 + wmargin, nwf, orbital)
# print("w=", wfplotmax, max(ywf))
kwf = wfplotmax / max(ywf)
ywf = [E + ywf[i] * kwf for i in range(len(ywf))]
WFlist_.append(ywf)
Siib, Siiw, Sii, Kb, Kw, K, Ub, Uw, U, Eb, Ew, E = cal_E_analytical(orbital)
print(" Sii = {:8.4f}(b) + {:8.4f}(w) = {:8.4f}".format(Siib, Siiw, Sii))
print(" K = {:8.4f}(b) + {:8.4f}(w) = {:8.4f} eV".format(Kb, Kw, K))
print(" U = {:8.4f}(b) + {:8.4f}(w) = {:8.4f} eV".format(Ub, Uw, U))
print(" E = {:8.4f}(b) + {:8.4f}(w) = {:8.4f} eV".format(Eb, Ew, E))
Sii, K, U, E = cal_E(orbital, xminEinteg, xmaxEinteg, nxEinteg)
print(" Sii={:8.4f}, K, U, E = ".format(Sii), K, U, E, "eV")
print("")
atoms[iatom]["orbitals"] = orbitals
atoms[iatom]["xwf"] = xwf
atoms[iatom]["wfs"] = WFlist_
# atom_orbitals.append(orbitals)
xWFlist.append(xwf)
yWFlist.append(WFlist_)
fig, ax = plt.subplots(1, 1)
ax1 = ax #ax[0]
# ax2をax1に関連させる
# ax2 = ax1.twinx()
ax1.plot(xpot, ypot, linewidth = 0.5) #, label = '$V$($x$)')
# Plot energy levels for infinite well potential
for i in range(nsites):
name, iatom, atom, x, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
for iL in range(1, 100):
E = quantization_level_inifinite(Vbase, meff, wwidth, iL)
if E > Vvac:
break
ax1.plot([xwmin, xwmax], [E, E], linestyle = 'dotted', color = 'r', linewidth = 0.5)
# Plot energy levels and wave functions for finite well potential
for i in range(nsites):
name, iatom, atom, x, wwidth, xwmin, xwmax, Vbase, meff = siteinf(i, crystal)
# orbitals = atom_orbitals[i]
orbitals = atoms[iatom]["orbitals"]
norbitals = len(orbitals)
for ie in range(norbitals):
orbital = orbitals[ie]
E = orbital["E"]
ax1.plot([xwmin, xwmax], [E, E], linestyle = 'dashed', color = 'b', linewidth = 1.0)
ax1.plot(xWFlist[i], yWFlist[i][ie], linestyle = '-', color = 'g', linewidth = 0.5)
ax1.set_xlabel("x ($\AA$)", fontsize = fontsize)
ax1.set_ylabel("$V$($x$), $\phi$($x$)", fontsize = fontsize)
ax1.legend(fontsize = legend_fontsize)
ax1.tick_params(labelsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
print("")
print("Press ENTER to exit>>", end = '')
input()
terminate()
[ドキュメント]
def main():
"""
プログラムのメインエントリポイントです。
グローバル変数`mode`の値に応じて、`view_basis()`, `band()`, `wf()`のいずれかの関数を呼び出します。
無効なモードが指定された場合はエラーメッセージを表示して終了します。
:returns: None
"""
global mode
if mode == 'basis':
view_basis()
elif mode == 'band':
band()
elif mode == 'wf':
wf()
else:
terminate("Error: Invalid mode [{}]".format(mode))
if __name__ == "__main__":
main()