anharmonic.py

プログラムの動作

このプログラムは、1次元単原子鎖（最近接相互作用モデル）における非調和フォノンの振動数が、温度変化によってどのようにシフトするかを計算するシミュレーションスクリプトです。

主な機能として、指定された波数パラメータ \(q\) に対して、指定された温度範囲で以下の値を計算します。

調和近似に基づく基本のフォノン振動数
三次の非調和性（K3）に由来する振動数シフト（Bubble図形に相当）
四次の非調和性（K4）に由来する振動数シフト（Tadpole図形に相当）

これらを摂動論的に足し合わせることで、温度に依存するフォノン振動数の変化を評価し、熱膨張やフォノン間相互作用などによる物理的な影響（非調和効果）を解析するという課題を解決します。

原理

1次元鎖の相対変位を \(\Delta_n = u_{n+1} - u_n\) としたとき、ポテンシャルエネルギー \(V\) は以下の式でモデル化されます。

\[V = \sum_n \left[ \frac{K2}{2} \Delta_n^2 + \frac{K3}{3!} \Delta_n^3 + \frac{K4}{4!} \Delta_n^4 \right]\]

ここで、\(K2\) は調和力定数、\(K3\) および \(K4\) はそれぞれ三次・四次の非調和定数です。調和近似における分散関係（振動数 \(\omega_q\)）は以下の式で与えられます。

\[\omega_q = \sqrt{ \frac{4 K2}{M} \sin^2\left(\frac{q a}{2}\right) }\]

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

\[\omega_q(T) \approx \omega_q + \Delta\omega_q^{(4)}(T) + \Delta\omega_q^{(3)}(T)\]

各シフト量は以下の原理に基づいて計算されます。

四次の寄与（Tadpoleシフト） \(\Delta\omega_q^{(4)}\): 熱平均された変位の二乗 \(\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)\]

三次の寄与（Bubbleシフト） \(\Delta\omega_q^{(3)}\): 波数空間 \(k\) に対する和（主値積分的な計算）として評価されます。計算の際、分母がゼロになる特異点での発散を防ぐため、微小な正則化パラメータ \(\eta\) を導入し、以下の式を利用して実部を抽出します。

\[\operatorname{Re}\left[ \frac{1}{x + i\eta} \right] = \frac{x}{x^2 + \eta^2}\]

これらの計算では、デフォルトで量子力学的なBose-Einstein分布関数が用いられますが、オプションによって古典極限の分布に切り替えることも可能です。

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

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

pip install numpy

必要な入力ファイル

本プログラムは、外部からの入力ファイルを必要としません。計算に使用される力定数（\(K2, K3, K4\)）は、ソースコード内に定義されている関数（k2_func, k3_func, k4_func）の戻り値として直接記述されています。特定の力定数やその温度依存性を設定したい場合は、これらの関数内の数値を直接書き換えて使用することが想定されています。

生成される出力ファイル

計算結果はCSVファイルとして保存されます。デフォルトの保存ファイル名は omega_vs_T.csv です。

出力されるCSVファイルには以下の内容（列）が含まれます。

T: 温度
omega0: 調和近似におけるフォノン振動数
domega3_bubble: 三次の非調和性による振動数シフト量
domega4_tadpole: 四次の非調和性による振動数シフト量
omegaT: 最終的な非調和フォノン振動数（omega0 + domega3_bubble + domega4_tadpole）

デフォルトでは振動数の単位は rad/s となりますが、コマンドライン引数の指定によって THz に変換された値が出力されます。

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

基本的な実行コマンドは以下の通りです。

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

利用可能な主な引数は以下の通りです。

--a: 格子定数 [m]（デフォルト: 1.0）
--M: 質量 [kg]（デフォルト: 1.0）
--Nk: ブリルアンゾーン内で分割する波数 \(k\) のグリッド数（デフォルト: 4096）
--kcut: 赤外発散（IR発散）を避けるためのカットオフ波数。\(1/a\) 単位（デフォルト: 0.0）
--q: 評価対象となる波数 \(q\) を \(\pi/a\) に対する比率で指定（デフォルト: 0.5）
--Tmin: 計算する温度の最小値 [K]（デフォルト: 0.0）
--Tmax: 計算する温度の最大値 [K]（デフォルト: 300.0）
--nT: 指定した温度範囲内で計算する点数（デフォルト: 31）
--classical: 量子統計（Bose分布）の代わりに、古典極限（高温極限）の統計を使用するフラグ
--eta_rel: Bubbleシフト計算時の正則化パラメータ \(\eta\) を決める係数（\(\eta = \text{eta\_rel} \times \omega_0\)）（デフォルト: 1e-3）
--out: 出力するCSVファイル名（デフォルト: omega_vs_T.csv）
--thz: 出力される振動数の単位を rad/s から THz に変換するフラグ

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

python D:/git/sphinx/tkProg/source/tiny_simulations/anharmonic.py --q 1.0 --Tmin 10 --Tmax 500 --nT 50 --thz --out result.csv

実行結果の説明: このコマンドを実行すると、波数 \(q = 1.0 \times (\pi/a)\)（ブリルアンゾーン境界）におけるフォノン振動数の温度依存性を評価します。温度は 10 K から 500 K までの範囲を 50 等分した各点で計算が行われます。 --thz オプションが指定されているため、計算された振動数とシフト量はすべて THz 単位に変換されます。結果はカレントディレクトリの result.csv というファイルに保存され、処理完了時には標準出力にファイル名と書き込まれた列のヘッダー情報、および単位（THz）が表示されます。