(section_1_title)="1. ドキュメントタイトル"

pw1d.py プログラムドキュメント

(section_2_overview)="2. 概要" このドキュメントは、Pythonスクリプト pw1d.py の機能と使用法を解説します。

pw1d.py は、1次元の平面波基底を用いてバンド計算を行うスクリプトです。自由電子モデルに周期ポテンシャルを導入し、固体中の電子のエネルギーバンド構造や波動関数を計算・可視化します。FFTを利用してポテンシャルのフーリエ変換を計算し、平面波基底を用いて固有方程式を解きます。

関連リンク: コードのDocstringに :doc:`pw1d_usage``` と記述されていますが、このドキュメントセットには pw1d_usage`` というドキュメントが含まれていないため、参照エラーを避けるためにリンクは解除します。

(section_3_dependencies)="3. 依存ライブラリ" pw1d.py は以下のPythonライブラリに依存しています。

  • 標準ライブラリ:

    • sys: システム固有のパラメータと機能にアクセスするためのモジュール。

    • csv: CSV (Comma Separated Values) ファイルの読み書きをサポートするモジュール。

  • 非標準ライブラリ:

    • numpy (np): 数値計算のための基盤となるパッケージ。配列オブジェクトや線形代数ルーチンなどを提供します。

      • sqrt, exp, sin, cos, tan, pi: numpy から個別にインポートされる数学関数および定数。

      • numpy.linalg (LA): 線形代数ルーチンを含むモジュール。

    • matplotlib.pyplot (plt): グラフ描画のためのライブラリ。

これらのライブラリは、スクリプトを実行する前にインストールされている必要があります。

(section_4_configurations)="4. 設定とグローバル変数" pw1d.py は、スクリプトの冒頭で定義されるいくつかのグローバル変数によって動作を制御します。これらの一部はコマンドライン引数によって上書き可能です。

(subsection_4_1_physical_constants)="4.1. 物理定数" プログラム内で使用される物理定数です。

  • 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^2 N^-1 m^-2)

  • kB: ボルツマン定数 (1.380658e-23 JK^-1)

  • 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)。単位はeVとÅngströmに合わせて調整されています。

(subsection_4_2_global_configuration)="4.2. グローバル設定" プログラムの実行モードを制御します。

  • mode: 実行モード ('ft' | 'band' | 'wf')。デフォルトは 'ft'。

(subsection_4_3_crystal_definition)="4.3. 結晶定義" 結晶の特性を定義します。

  • a: 格子定数 (5.4064 Å)。デフォルト値はシリコンの格子定数です。

  • na: FFTのための分割数 (64)。2のべき乗である必要があります。

(subsection_4_4_potential)="4.4. ポテンシャル" ポテンシャルの種類とパラメータを定義します。

  • pottype: ポテンシャルの種類 ('rect' | 'gauss')。デフォルトは 'rect'。

  • bwidth: 障壁幅 (0.5 Å)。

  • bpot: 障壁高さ (10.0 eV)。

(subsection_4_5_band)="4.5. バンド計算設定" バンド計算のパラメータを定義します。

  • nG: G点数（基底関数の数） (3)。

  • kmin: k点の最小値 (-0.5)。単位は pi/a の係数。

  • kmax: k点の最大値 (0.5)。単位は pi/a の係数。

  • nk: k点の分割数 (21)。

  • Erange: プロットするエネルギー範囲 ([0.0, 10.0] eV)。

(subsection_4_6_wave_function)="4.6. 波動関数計算設定" 波動関数計算のパラメータを定義します。

  • xwmin: 波動関数プロットのx座標最小値 (0.0 Å)。

  • xwmax: 波動関数プロットのx座標最大値 (3.0 * a Å)。

  • nxw: 波動関数プロットのx点数 (101)。

  • kw: 波動関数を計算するk値 (0.0)。単位は pi/a の係数。

  • iLevel: プロットするエネルギー準位のインデックス (0)。

(subsection_4_7_figure_configuration)="4.7. 図の設定" Matplotlibによる図の表示設定です。

  • figsize: 図のサイズ ((6, 8))。

  • fontsize: フォントサイズ (12)。

  • legend_fontsize: 凡例のフォントサイズ (12)。

(section_5_command_line_arguments)="5. コマンドライン引数" プログラムはコマンドライン引数を受け取ることができ、これにより上記グローバル変数のデフォルト値を上書きします。引数はモードによって必要な数が異なります。

基本的な呼び出し形式は以下の通りです。

python pw1d.py (mode a na pottype bwidth bpot ...)

引数の順序と意味は以下の通りです。括弧内の引数はオプションであり、指定しない場合はデフォルト値が使用されます。

  • 全モード共通:

    1. mode (str): 実行モード ('ft', 'band', 'wf')。デフォルトは 'ft'。

    2. a (float): 格子定数 (Å)。デフォルトは 5.4064。

    3. na (int): FFTのための分割数。デフォルトは 64。

    4. pottype (str): ポテンシャルの種類 ('rect', 'gauss')。デフォルトは 'rect'。

    5. bwidth (float): 障壁幅 (Å)。デフォルトは 0.5。

    6. bpot (float): 障壁高さ (eV)。デフォルトは 10.0。

  • band モードの場合 (上記共通引数に続けて): 7. nG (int): G点数（基底関数の数）。デフォルトは 3。 8. kmin (float): k点の最小値 (pi/a 単位)。デフォルトは -0.5。 9. kmax (float): k点の最大値 (pi/a 単位)。デフォルトは 0.5。 10. nk (int): k点の分割数。デフォルトは 21。

  • wf モードの場合 (上記共通引数に続けて): 7. nG (int): G点数（基底関数の数）。デフォルトは 3。 8. kw (float): 波動関数を計算するk値 (pi/a 単位)。デフォルトは 0.0。 9. iLevel (int): プロットするエネルギー準位のインデックス。デフォルトは 0。 10. xwmin (float): 波動関数プロットのx座標最小値 (Å)。デフォルトは 0.0。 11. xwmax (float): 波動関数プロットのx座標最大値 (Å)。デフォルトは 3.0 * a。 12. nxw (int): 波動関数プロットのx点数。デフォルトは 101。

引数が不足している場合、または不正な場合は、デフォルト値が適用されるか、usage() メッセージと共にプログラムが終了します。

使用例:

python pw1d.py
python pw1d.py ft 5.4064 64 rect 0.5 10.0
python pw1d.py band 5.4064 64 rect 0.5 10.0 3 -0.5 0.5 21
python pw1d.py wf 5.4064 64 rect 0.5 10.0 3 0.0 0 0.0 16.2192 101

(section_6_function_descriptions)="6. 関数説明" (subsection_6_1_pfloat)="6.1. pfloat(str)"

  • 概要: 文字列を浮動小数点数に安全に変換します。

  • 詳細: 変換できない場合は None を返し、プログラムは終了しません。

  • パラメータ:

    • str (str): 変換する文字列。

  • 戻り値: float または None: 変換された浮動小数点数、または変換できなかった場合は None。

(subsection_6_2_pint)="6.2. pint(str)"

  • 概要: 文字列を整数に安全に変換します。

  • 詳細: 変換できない場合は None を返し、プログラムは終了しません。

  • パラメータ:

    • str (str): 変換する文字列。

  • 戻り値: int または None: 変換された整数、または変換できなかった場合は None。

(subsection_6_3_getarg)="6.3. getarg(position, defval=None)"

  • 概要: コマンドライン引数を安全に取得します。

  • 詳細: 指定されたインデックスの引数が存在しない場合は、デフォルト値を返します。

  • パラメータ:

    • position (int): 取得する引数の位置 (sys.argv のインデックス)。

    • defval (any, optional): 引数が存在しない場合に返すデフォルト値。デフォルトは None。

  • 戻り値: str または any: 取得した引数の文字列、またはデフォルト値。

(subsection_6_4_getfloatarg)="6.4. getfloatarg(position, defval=None)"

  • 概要: コマンドライン引数を浮動小数点数として安全に取得します。

  • 詳細: getarg と pfloat を組み合わせて、引数を浮動小数点数に変換して返します。

  • パラメータ:

    • position (int): 取得する引数の位置。

    • defval (any, optional): 引数が存在しない場合、または変換できなかった場合に返すデフォルト値。デフォルトは None。

  • 戻り値: float または any: 変換された浮動小数点数、またはデフォルト値。

(subsection_6_5_getintarg)="6.5. getintarg(position, defval=None)"

  • 概要: コマンドライン引数を整数として安全に取得します。

  • 詳細: getarg と pint を組み合わせて、引数を整数に変換して返します。

  • パラメータ:

    • position (int): 取得する引数の位置。

    • defval (any, optional): 引数が存在しない場合、または変換できなかった場合に返すデフォルト値。デフォルトは None。

  • 戻り値: int または any: 変換された整数、またはデフォルト値。

(subsection_6_6_usage)="6.6. usage()"

  • 概要: スクリプトの正しい使用方法を標準出力に表示します。

  • 詳細: コマンドライン引数の形式と、利用可能なモード (ft, band, wf) およびそれに対応するパラメータの例を示します。

  • パラメータ: なし。

  • 戻り値: なし。

(subsection_6_7_terminate)="6.7. terminate(message=None)"

  • 概要: メッセージを表示し、使用方法を表示してプログラムを終了します。

  • 詳細: エラーが発生した場合などに呼び出され、ユーザーに情報を提供し、スクリプトの実行を停止します。

  • パラメータ:

    • message (str, optional): 終了前に表示するエラーメッセージ。デフォルトは None。

  • 戻り値: なし (プログラムが終了するため)。

(subsection_6_8_reduce)="6.8. reduce(x, x0)"

  • 概要: x を x0 の周期で還元します。

  • 詳細: x を [0, x0) の範囲にマッピングします。負の値の場合も正しく処理されます。

  • パラメータ:

    • x (float): 還元する値。

    • x0 (float): 周期。

  • 戻り値: float: 還元された値。

(subsection_6_9_pot)="6.9. pot(x)"

  • 概要: 指定された位置でのポテンシャルエネルギーを計算します。

  • 詳細: グローバル変数 pottype によって矩形障壁 (rect) またはガウス型ポテンシャル (gauss) のいずれかを計算します。周期境界条件を考慮して reduce 関数を使用します。

  • パラメータ:

    • x (float): ポテンシャルを評価する空間座標。

  • 戻り値: float: x におけるポテンシャルエネルギー (eV)。

(subsection_6_10_build_potential)="6.10. build_potential(xmin, xstep, n)"

  • 概要: 空間グリッド上のポテンシャル分布を構築します。

  • 詳細: pot 関数を使用して、指定された範囲とステップでポテンシャル値を計算し、x座標とポテンシャル値の配列を返します。

  • パラメータ:

    • xmin (float): ポテンシャル計算を開始するx座標。

    • xstep (float): x座標のステップサイズ。

    • n (int): 計算する点の数。

  • 戻り値: tuple[numpy.ndarray, numpy.ndarray]:

    • x座標の配列。

    • 対応するポテンシャル値の配列。

(subsection_6_11_find_iG)="6.11. find_iG(dij, Glist)"

  • 概要: 指定された逆格子座標 dij に対応する Glist のインデックスを検索します。

  • 詳細: Glist 内の要素と dij を比較し、一致する要素のインデックスを返します。

  • パラメータ:

    • dij (int): 検索する逆格子座標。

    • Glist (list[int] or numpy.ndarray[int]): 逆格子座標のリスト。

  • 戻り値: int または None: dij が見つかった場合のインデックス、見つからなかった場合は None。

(subsection_6_12_find_Vft)="6.12. find_Vft(dij, Glist, Vft)"

  • 概要: 指定された逆格子座標 dij に対応するフーリエ変換されたポテンシャル値 Vft を検索します。

  • 詳細: Glist を使用して dij のインデックスを見つけ、そのインデックスに対応する Vft の値を返します。見つからなかった場合は 0.0 を返します。

  • パラメータ:

    • dij (int): 検索する逆格子座標。

    • Glist (list[int] or numpy.ndarray[int]): 逆格子座標のリスト。

    • Vft (list[complex] or numpy.ndarray[complex]): フーリエ変換されたポテンシャル値のリスト。

  • 戻り値: complex または float: dij に対応する Vft の値、または 0.0。

(subsection_6_13_cal_fft)="6.13. cal_fft(na, a, ypot)"

  • 概要: 実空間ポテンシャルのフーリエ変換を計算し、物理的な逆格子空間の表現に変換します。

  • 詳細: numpy のFFT関数を使用し、結果をG<0の成分が負のG値に対応するようにシフトします。

  • パラメータ:

    • na (int): FFTに用いるサンプリング点の数。

    • a (float): 格子定数（単位長さ）。

    • ypot (numpy.ndarray[float]): 実空間のポテンシャル値の配列。

  • 戻り値: tuple:

    • xft0 (list[int]): FFTの生の結果のx軸インデックス (0から na-1)。

    • yft0 (numpy.ndarray[complex]): FFTの生の結果のフーリエ変換されたポテンシャル。

    • xft (numpy.ndarray[float]): 物理的に並べ替えられた逆格子座標 (G値)。

    • yft (numpy.ndarray[complex]): 物理的に並べ替えられたフーリエ変換されたポテンシャル。

    • iGlist (list[int]): G値に対応する逆格子点の整数インデックス。

    • nahalf (int): FFT点の数の半分。

    • xftstep (float): 逆格子空間のステップサイズ。

(subsection_6_14_extract_basis)="6.14. extract_basis(yft0, nG)"

  • 概要: フーリエ変換されたポテンシャルから、指定された数のG点基底（平面波基底）を抽出します。

  • 詳細: nG に基づき、中心 (G=0) から対称的にG点を抽出します。nG=2 の場合はG=0とG=1を、それ以外の場合はG=0を中心に正負のG点を抽出します。抽出された Glist の要素の順序は、昇順にはソートされず、[0, 1, ..., nGmax, -1, ..., -nGmax] のような順序になります。

  • パラメータ:

    • yft0 (numpy.ndarray[complex]): FFTの生の結果のフーリエ変換されたポテンシャル。

    • nG (int): 抽出するG点の総数。

  • 戻り値: tuple:

    • Glist (numpy.ndarray[int]): 抽出されたG点の整数インデックスのリスト。

    • Vftlist (numpy.ndarray[complex]): 抽出されたG点に対応するフーリエ変換されたポテンシャル値のリスト。

    • nGmax (int): 抽出されたG点の最大絶対値 (int(nG / 2))。

(subsection_6_15_free_KE)="6.15. free_KE(k, iG)"

  • 概要: 自由電子の運動エネルギーを計算します。

  • 詳細: ブロッホ波の波数 k と逆格子ベクトル iG を用いて、運動エネルギーを求めます。

  • パラメータ:

    • k (float): ブロッホ波の波数。単位は pi/a の係数。

    • iG (int): 逆格子点の整数インデックス。

  • 戻り値: float: 計算された運動エネルギー (eV)。

(subsection_6_16_solve_pw)="6.16. solve_pw(Glist, Vftlist, nGmax, nG, k, IsPrint=0)"

  • 概要: 平面波基底を用いて固有方程式を解き、エネルギー固有値と平面波係数を計算します。

  • 詳細: フォック行列（ハミルトニアン行列）を構築し、Numpy の線形代数ソルバーを用いて固有値問題（Schrödinger方程式）を解きます。ハミルトニアンはエルミート行列となるように構築されます。

  • パラメータ:

    • Glist (numpy.ndarray[int]): 考慮する逆格子点 (G値) のリスト。

    • Vftlist (numpy.ndarray[complex]): 各G点に対応するフーリエ変換されたポテンシャル。

    • nGmax (int): Glist における最大の正のG値。

    • nG (int): 基底関数の総数。

    • k (float): 計算を行うブロッホ波の波数。

    • IsPrint (int, optional): デバッグ情報の出力レベル。0 はなし、1 は主要な値、2 は詳細な値を出力。デフォルトは 0。

  • 戻り値: tuple[numpy.ndarray[float], numpy.ndarray[complex]]:

    • ei (numpy.ndarray[float]): 計算されたエネルギー固有値の配列。

    • ci (numpy.ndarray[complex]): 各エネルギー固有値に対応する平面波係数（固有ベクトル）。

(subsection_6_17_cal_wf)="6.17. cal_wf(xwmin, xwmax, nxw, kw, ci, Glist)"

  • 概要: 平面波係数から実空間の波動関数を計算します。

  • 詳細: 指定された空間範囲とグリッドで、平面波の重ね合わせとして波動関数とその確率密度を計算します。

  • パラメータ:

    • xwmin (float): 波動関数を計算するx座標の最小値。

    • xwmax (float): 波動関数を計算するx座標の最大値。

    • nxw (int): 計算するx点の数。

    • kw (float): ブロッホ波の波数。

    • ci (list[complex] or numpy.ndarray[complex]): 各G点に対応する平面波の係数。

    • Glist (list[int] or numpy.ndarray[int]): 考慮する逆格子点 (G値) のリスト。

  • 戻り値: tuple[numpy.ndarray[complex], numpy.ndarray[complex], complex]:

    • xwf (numpy.ndarray[complex]): 波動関数を評価したx座標の配列。

    • ywf (numpy.ndarray[complex]): 計算された波動関数の複素数値の配列。

    • charge_scalar (complex): 波動関数の確率密度 |Ψ|^2 の値。注意: この関数は charge として単一のスカラー値（ループの最後の x における確率密度）を返します。呼び出し元である wf() 関数では、この戻り値は使用されず、 ywf から確率密度が配列として再計算されます。

(subsection_6_18_wf)="6.18. wf()"

  • 概要: 波動関数を計算し、ポテンシャルと合わせてプロットします。

  • 詳細: 設定されたパラメータに基づき、ポテンシャル、フーリエ変換されたポテンシャル、固有値、選択されたエネルギー準位の波動関数を計算し、Matplotlib で可視化します。計算結果は標準出力にも表示されます。

  • パラメータ: なし。

  • 戻り値: なし。

(subsection_6_19_band)="6.19. band()"

  • 概要: エネルギーバンド構造を計算し、自由電子バンドと共にプロットします。

  • 詳細: 異なる波数 k に対して固有方程式を解き、エネルギー固有値を計算します。計算されたバンド構造は、自由電子のエネルギーバンドと比較して Matplotlib でプロットされます。nG が 4 以上の場合、自動的に奇数に変換されます (例: 4 -> 3, 5 -> 5, 6 -> 5)。nG の値が 0 以下または 64 より大きい場合はエラーで終了します。

  • パラメータ: なし。

  • 戻り値: なし。

(subsection_6_20_ft)="6.20. ft()"

  • 概要: ポテンシャルのフーリエ変換を計算し、実空間ポテンシャルとフーリエ変換されたポテンシャルをプロットします。

  • 詳細: 設定されたパラメータに基づいて実空間ポテンシャルを構築し、そのフーリエ変換を実行します。結果は Matplotlib を用いて、実空間と逆格子空間の両方で可視化されます。

  • パラメータ: なし。

  • 戻り値: なし。

(subsection_6_21_main)="6.21. main()"

  • 概要: スクリプトの主要な実行ロジック。

  • 詳細: グローバル変数 mode に基づいて、ft(), band(), wf() のいずれかの関数を呼び出します。無効なモードが指定された場合はエラーメッセージと共に終了します。

  • パラメータ: なし。

  • 戻り値: なし。

(section_7_execution)="7. 実行方法" pw1d.py はPythonインタプリタを使用して実行されます。コマンドライン引数により動作モードやパラメータを制御できます。

python pw1d.py [オプション引数...]

オプション引数の例:

  • デフォルトモード (ft):

    python pw1d.py

  • フーリエ変換モード (ft) で特定のパラメータを指定:

    python pw1d.py ft 5.0 128 rect 0.3 5.0

  • バンド計算モード (band) で特定のパラメータを指定:

    python pw1d.py band 5.4064 64 rect 0.5 10.0 5 -0.3 0.3 31

  • 波動関数モード (wf) で特定のパラメータを指定:

    python pw1d.py wf 5.4064 64 rect 0.5 10.0 5 0.1 1 0.0 16.2192 201

プログラムの実行後、グラフが表示され、ユーザーが ENTER キーを押すまで一時停止します。

(section_8_input_output)="8. 入出力" (subsection_8_1_input)="8.1. 入力"

  • コマンドライン引数: プログラムの動作モード (mode) や各種物理的・計算的パラメータは、コマンドライン引数として入力されます。

  • ユーザー入力: グラフ表示後、プログラムはユーザーが ENTER キーを押すのを待ちます。これは input() 関数によって処理されます。

pw1d.py はファイルからの入力を直接サポートしていません。すべての設定はスクリプト内のグローバル変数またはコマンドライン引数を通じて行われます。

(subsection_8_2_output)="8.2. 出力"

  • 標準出力 (stdout):

    • プログラムの起動時、設定されたパラメータが表示されます。

    • FFTの結果、G基底の抽出、固有方程式の解（エネルギー固有値や平面波係数）などの計算の中間結果や詳細情報が表示されます。

    • エラーメッセージや使用方法 (usage()) が表示されることがあります。

  • グラフ表示 (matplotlib.pyplot):

    • ft モード: 実空間ポテンシャルと、そのフーリエ変換されたポテンシャル（生の結果と物理的に並べ替えられた結果）がグラフとして表示されます。

    • band モード: 計算されたエネルギーバンド構造が、自由電子のエネルギーバンドと比較してプロットされます。

    • wf モード: 選択されたエネルギー準位の波動関数（実部、虚部）と確率密度、および実空間ポテンシャルがプロットされます。

    • すべてのグラフは plt.pause(0.1) の後、ユーザーが ENTER を押すまで表示が維持されます。

pw1d.py は計算結果をファイルに保存する機能は持っていません。

(section_9_notes)="9. 注意事項"

  • cal_wf() の charge 戻り値:

    • cal_wf() 関数の charge 戻り値は、Docstringでは「波動関数の確率密度」と記述されていますが、実際にはループの最後の x 座標における単一のスカラー値 (f * f.conjugate()) が返されます。

    • この関数の呼び出し元である wf() 関数では、この戻り値は利用されず、 ywf 配列から charge = [(f0 * f0.conjugate()).real for f0 in ywf] として確率密度が配列として再計算されプロットされます。

  • band() モードにおける nG の調整:

    • band() 関数内では、nG (基底関数の数) が 4 以上の場合、nG = int(nG / 2) * 2 + 1 というロジックで奇数に調整されます。これにより、G=0 を中心とした対称的なG点のセットが維持されやすくなります (例: nG=4 -> 3, nG=6 -> 5)。

    • nG が 0 以下または 64 より大きい場合はエラーで終了します。

  • グラフの凡例:

    • band() 関数および ft() 関数内のいくつかの plot() 呼び出しでは、label 引数が指定されていないため、ax1.legend() や ax1.legend(fontsize = legend_fontsize) が呼び出されても凡例は表示されません。

  • ft() 関数の cal_fft 戻り値の変数名:

    • ft() 関数内で cal_fft の戻り値を受け取る際に、iGlist を意味する戻り値が iG という変数名で受け取られています。しかし、iG はその後のコードではリストとして適切に扱われています。

  • プロットのX軸範囲:

    • band() 関数において、プロットのX軸範囲は ax1.set_xlim([-0.5, 0.5]) とハードコードされており、コマンドライン引数の kmin, kmax の値に連動していません。