D:/git/sphinx/tkProg/source/tiny_simulations/anharmonic3D.py

プログラムの動作

本プログラムは、1次元単原子連鎖モデルにおける非調和フォノンの振動数シフトの温度依存性を計算するシミュレーションスクリプトです。 隣接原子間の相互作用に3次および4次の非調和項（それぞれ \(K_3\), \(K_4\)）を導入し、摂動論を用いて特定の波数 \(q\) におけるフォノンの自己エネルギーを計算します。 これにより、温度上昇に伴うフォノンの熱的な軟化（または硬化）を定量的に評価し、非調和性が格子力学に与える影響を解析するという課題を解決します。 力の定数はソースコード内の関数を編集することで、ユーザーが任意の関数として定義することが可能です。

原理

本モデルでは、隣接する原子間の相対変位 \(\Delta_n = u_{n+1} - u_n\) に対して、以下のポテンシャルエネルギーを仮定しています。

\[ V = \sum_n \left[ \frac{K_2}{2} \Delta_n^2 + \frac{K_3}{3!} \Delta_n^3 + \frac{K_4}{4!} \Delta_n^4 \right] \]

調和近似におけるフォノンの分散関係は、格子定数を \(a\)、質量を \(M\) として以下の式で与えられます。

\[ \omega_q = \sqrt{\frac{4 K_2}{M} \sin^2\left(\frac{qa}{2}\right)} \]

非調和性による振動数シフト \(\Delta\omega_q(T)\) は、自己エネルギー \(\Sigma(q, \omega_q)\) の実部を用いて次のように近似されます。

\[ \Delta\omega_q(T) \approx \frac{\text{Re} \Sigma(q, \omega_q)}{2 \omega_q} \]

本プログラムでは、以下の2つの寄与を計算します。

  1. 4次非調和項（Tadpoleシフト）: 4次定数 \(K_4\) に由来するシフトで、他の波数 \(k\) の熱平均 \(\langle |u_k|^2 \rangle\) を用いて計算されます。変位の自乗平均は、量子論的な揺らぎと熱的な揺らぎを含み、次のように表されます。

\[ \langle |u_k|^2 \rangle = \frac{\hbar}{2 M \omega_k} \coth\left(\frac{\hbar \omega_k}{2 k_B T}\right) \]

  1. 3次非調和項（Bubbleシフト）: 3次定数 \(K_3\) に由来するシフトで、運動量保存則を満たす2つのフォノンへの崩壊と再結合の過程を表します。エネルギー分母に微小な正則化パラメータ \(\eta\) を導入し、主値積分を用いて計算されます。ボース分布関数 \(n_k\) を用いて、各モード間の散乱確率を重み付けしています。

最終的な温度 \(T\) での振動数は、\(\omega_q(T) \approx \omega_q + \Delta\omega_q^{(3)}(T) + \Delta\omega_q^{(4)}(T)\) として求められます。

必要な非標準ライブラリとインストール方法

本プログラムは数値計算ライブラリである numpy を使用しています。実行前に以下のコマンドを用いてインストールしてください。

\[ \text{pip install numpy} \]

その他の argparse, math, dataclasses, typing などはPython標準ライブラリであるため、追加のインストールは不要です。

必要な入力ファイル

本プログラムの実行にあたり、外部からの入力ファイルは必要ありません。 モデルにおける各力の定数（\(K_2\), \(K_3\), \(K_4\)）は、ソースコード内で定義されている関数 k2_func, k3_func, k4_func によって決定されます。 条件を変更したい場合は、コード中の対象関数（3重引用符で説明が記述されている部分の直下）の戻り値を直接編集してください。

生成される出力ファイル

プログラムを実行すると、指定したファイル名（デフォルトは omega_vs_T.csv）で計算結果がCSV形式として保存されます。

出力ファイルには以下の5列のデータが記録されます。

  • T: 温度 [K]

  • omega0: 調和振動数 \(\omega_q\) （温度非依存）

  • domega3_bubble: 3次非調和項による振動数シフト \(\Delta\omega_q^{(3)}\)

  • domega4_tadpole: 4次非調和項による振動数シフト \(\Delta\omega_q^{(4)}\)

  • omegaT: 最終的な非調和振動数 \(\omega_q(T)\)

出力される振動数の単位は、デフォルトでは rad/s ですが、オプションを指定することで THz に変更可能です。

コマンドラインでの使用例 (Usage)

本プログラムはコマンドライン引数を利用して、計算のパラメータを柔軟に変更できます。基本的な実行コマンドは以下の通りです。

python D:/git/sphinx/tkProg/source/tiny_simulations/anharmonic3D.py [オプション]

利用可能な主なオプション引数は以下の通りです。

  • --a : 格子定数 [m]（デフォルト: 1.0）

  • --M : 質量 [kg]（デフォルト: 1.0）

  • --Nk : ブリルアンゾーン内のk点数（デフォルト: 4096）

  • --kcut : 赤外発散を回避するために除外する波数の閾値 [1/a]（デフォルト: 0.0）

  • --q : 計算対象となる波数を \(\pi/a\) の割合で指定（デフォルト: 0.5）

  • --Tmin : 最小温度 [K]（デフォルト: 0.0）

  • --Tmax : 最大温度 [K]（デフォルト: 300.0）

  • --nT : 計算する温度のステップ数（デフォルト: 31）

  • --classical : 量子統計の代わりに古典統計（高温極限）を使用する場合に指定

  • --eta_rel : バブルシフトの主値積分における正則化パラメータ \(\eta\) の係数（デフォルト: 1e-3）

  • --out : 出力するCSVのファイル名（デフォルト: omega_vs_T.csv）

  • --thz : 振動数の出力単位を rad/s ではなく THz (cycles/s) にする場合に指定

コマンドラインでの具体的な使用例

温度範囲を 0 K から 500 K まで 51 ステップで計算し、出力単位を THz にして result_phonon.csv というファイルに保存する実行例です。

python D:/git/sphinx/tkProg/source/tiny_simulations/anharmonic3D.py --Tmax 500.0 --nT 51 --thz --out result_phonon.csv

実行結果の説明 コマンドを実行すると、ターミナルに以下のようなメッセージが表示され、計算が正常に終了したことが通知されます。

Wrote: result_phonon.csv
Columns: T,omega0,domega3_bubble,domega4_tadpole,omegaT
Frequencies are in THz.

カレントディレクトリに生成された result_phonon.csv を開くと、0 K から 500 K までの温度依存性が記述されており、各温度における3次非調和による負のシフトと、4次非調和による正のシフトの競合関係や、温度上昇に伴う非調和フォノン振動数の変化を THz 単位で確認することができます。