2025-05-25

数演算子のフーリエ変換と波数空間の生成・消滅演算子のペアとの違い

量子多体系の議論で、「（電子）数演算子のフーリエ変換」が登場するが、パッと見で「波数空間での生成・消滅演算子のペア」と同じなのか違うのかがわかりにくい。
ここでメモしておく。

格子上での（電子）数演算子は、位置空間（サイト表示）で次のように表される
$\displaystyle n_j = c_j^\dagger c_j$

これのフーリエ変換は
$\displaystyle \sqrt{N} \rho_q = \sum_j e^{-i q R_j} n_j = \sum_j e^{-i q R_j} c_j^\dagger c_j = \sum_k c_{k+q}^\dagger c_k$
となり、 $q=0$ のときに「生成・消滅演算子のペア」と同一になる。

途中の式変形は以下の通り。
$\displaystyle \sum_j e^{-i q R_j} c_j^\dagger c_j = \frac{1}{N} \sum_{j,k,k'} e^{-i q R_j} e^{i k' R_j} c_{k'}^\dagger e^{-i k R_j} c_{k} \\ \displaystyle \qquad = \sum_{k,k'} \delta_{k' - k, q} c_{k'}^\dagger c_{k} = \sum_k c_{k+q}^\dagger c_k$

$q = 0$ ：全電子数
$q = \pi$ ：電荷密度波（CDW）

応用例：電子–格子相互作用

電子–格子相互作用：
$\displaystyle H_{e\text{-}ph} \sim \sum_q g_q (a_q + a_{-q}^\dagger) \rho_q = \sum_{k,q} g_q (a_q + a_{-q}^\dagger) c_{k+q}^\dagger c_k$
$\rho_q$ が「電子が格子変形と相互作用して運動量を変える過程」を記述していることが分かる。

2024-08-12

動力学行列の正確な導出

固体物理

以下の代表的な固体物理の教科書を開いても、原子の変位をFourier変換ではなく、単一の平面波で天下りで導出しているため、一般性・正確性に欠ける。

キッテル
グロッソ、パラビチニ
イバッハ
パインズ

ここでは、キチンと原子の変位をFourier変換して動力学行列を導出する。

原子の変位を次のようにフーリエ変換する。
$\displaystyle \mathbf{u}_{\alpha}(\mathbf{R}_l, t) = \frac{1}{\sqrt{m_{\alpha}}} \frac{1}{N} \sum_{\mathbf{q}} \int_{-\infty}^{\infty} d \omega \, \mathbf{u}_{\alpha}(\mathbf{q}, \omega) e^{i (\mathbf{q} \cdot \mathbf{R}_l - \omega t)}$
$\frac{1}{\sqrt{m_{\alpha}}}$ は後の便宜のためであり、Fourier変換そのものとは関係が無い。

運動方程式は次のように表せる。
$\displaystyle m_{\alpha} \frac{d^2 \mathbf{u}_{\alpha}(\mathbf{R}_l, t)}{dt^2} = -\sum_{\beta, \mathbf{R}_{l'}} \Phi_{\alpha\beta}(\mathbf{R}_l - \mathbf{R}_{l'}) \mathbf{u}_{\beta}(\mathbf{R}_{l'}, t)$

ここで、 $\Phi$ と $\mathbf{r}$ の定義は以下の通り。
$\displaystyle \Phi_{\alpha\beta}(\mathbf{R}_l - \mathbf{R}_{l'}) = \left(\frac{ \partial V( \left\{ \mathbf{r} \right\} - \left\{ \mathbf{r}' \right\} ) }{ \partial \mathbf{r}_{\alpha l} \partial \mathbf{r}_{\beta l'} } \right)_{\{\mathbf{u}\}=0}$

$\displaystyle \mathbf{r}_{\alpha l}(t) = \mathbf{r}_{\alpha} + \mathbf{R}_l + \mathbf{u}_{\alpha}(\mathbf{R}_{l}, t)$

$V( \left\{ \mathbf{r} \right\} - \left\{ \mathbf{r}' \right\} )$ はポテンシャルが各粒子の相対位置にのみ依存することを示す。

ポイントとなるのは、結合定数 $\Phi$ のフーリエ変換をすることである。
$\displaystyle \Phi_{\alpha \beta}(\mathbf{R}_l - \mathbf{R}_{l'} ) = \frac{1}{N} \sum_{\mathbf{q}} \Phi_{\alpha \beta}(\mathbf{q}) e^{i \mathbf{q} \cdot (\mathbf{R}_l - \mathbf{R}_{l'}) }$

まず、運動方程式の左辺をフーリエ変換する。
$\displaystyle \frac{d^2 \mathbf{u}_{\alpha}(\mathbf{R}_l, t)}{dt^2} = \frac{1}{\sqrt{m_{\alpha}}} \frac{1}{N} \sum_{\mathbf{q}} \int_{-\infty}^{\infty} d\omega \, (-\omega^2) \mathbf{u}_{\alpha}(\mathbf{q}, \omega) e^{i (\mathbf{q} \cdot \mathbf{R}_l - \omega t)}$

$\displaystyle m_{\alpha} \frac{d^2 \mathbf{u}_{\alpha}(\mathbf{R}_l, t)}{dt^2} = \frac{ \sqrt{ m_{\alpha} } }{ N } \sum_{\mathbf{q}} \int_{-\infty}^{\infty} d\omega \, (-\omega^2) \mathbf{u}_{\alpha}(\mathbf{q}, \omega) e^{i (\mathbf{q} \cdot \mathbf{R}_l - \omega t)}$

次に運動方程式の右辺をフーリエ変換する。
$\displaystyle -\sum_{\beta, \mathbf{R}_{l'}} \Phi_{\alpha\beta}(\mathbf{R}_l - \mathbf{R}_{l'}) \mathbf{u}_{\beta}(\mathbf{R}_{l'}, t) \\ \displaystyle = - \frac{1}{N^2} \sum_{\beta, \mathbf{R}_{l'}} \sum_{\mathbf{q'}} \Phi_{\alpha\beta}(\mathbf{q'}) e^{i \mathbf{q'} \cdot ( \mathbf{R}_{l} - \mathbf{R}_{l'}) } \sum_{\mathbf{q}} \int_{-\infty}^{\infty} d\omega \, \frac{1}{\sqrt{m_{\beta}}} \mathbf{u}_{\beta}(\mathbf{q}, \omega) e^{i (\mathbf{q} \cdot \mathbf{R}_{l'} - \omega t)} \\ \displaystyle = - \frac{1}{N^2} \sum_{\beta} \sum_{\mathbf{q},\mathbf{q'}} \frac{1}{\sqrt{m_{\beta}}} \Phi_{\alpha\beta}(\mathbf{q'}) \int_{-\infty}^{\infty} d\omega \, \mathbf{u}_{\beta}(\mathbf{q}, \omega) e^{i \mathbf{q'} \cdot \mathbf{R}_{l} } e^{-i \omega t} \sum_{\mathbf{R}_{l'}} e^{i (\mathbf{q} - \mathbf{q'}) \cdot \mathbf{R}_{l'} } \\ \displaystyle = - \frac{1}{N} \sum_{\beta} \sum_{\mathbf{q},\mathbf{q'}} \frac{1}{\sqrt{m_{\beta}}} \Phi_{\alpha\beta}(\mathbf{q'}) \int_{-\infty}^{\infty} d\omega \, \mathbf{u}_{\beta}(\mathbf{q}, \omega) e^{i \mathbf{q'} \cdot \mathbf{R}_{l} } e^{-i \omega t} \delta_{\mathbf{q}, \mathbf{q'}} \\ \displaystyle = - \frac{1}{N} \sum_{\beta} \sum_{\mathbf{q}} \frac{1}{\sqrt{m_{\beta}}} \Phi_{\alpha\beta}(\mathbf{q}) \int_{-\infty}^{\infty} d\omega \, \mathbf{u}_{\beta}(\mathbf{q}, \omega) e^{i \mathbf{q} \cdot \mathbf{R}_{l} } e^{-i \omega t}$

ここで、以下の関係を使った。
$\displaystyle \sum_{\mathbf{R}_l} e^{-i \mathbf{q'} \cdot \mathbf{R}_l} e^{i \mathbf{q} \cdot \mathbf{R}_l} = N \delta_{\mathbf{q'}, \mathbf{q}}$

両辺の結果をまとめると、
$\displaystyle \sum_{\mathbf{q}} \int_{-\infty}^{\infty} d\omega \, \left( -\omega^2 \mathbf{u}_{\alpha}(\mathbf{q}, \omega) \right) e^{i (\mathbf{q} \cdot \mathbf{R}_l - \omega t)} \\ \displaystyle = \sum_{\mathbf{q}} \int_{-\infty}^{\infty} d\omega \, \left( -\sum_{\beta} \frac{1}{\sqrt{m_{\alpha} m_{\beta}}} \Phi_{\alpha\beta}(\mathbf{q}) \mathbf{u}_{\beta}(\mathbf{q}, \omega) \right) e^{i (\mathbf{q} \cdot \mathbf{R}_{l} -i \omega t) }$

ここまでくれば、両辺に平面波 $e^{-i ( \mathbf{q'} \cdot \mathbf{R}_l - \omega' t) }$ を掛けて $t$ の積分およびの $\mathbf{R}$ の総和を取ることで、 $q', \omega'$ の方程式に変換できる。（その後で、 $(q', \omega') \rightarrow (q, \omega)$ と書き直して、見栄えを良くする）

冒頭に言及した書籍では、変位 $u$ をFourier変換するのではなく、単一の平面波成分のみを持つとして式変形を進めている。
この記事でも採用している、ポテンシャルの2階微分まで拾う調和近似の範囲であれば、並進対称性を持つ条件で単一の $q$ のみが現れる方程式が得られる。
それはすなわち、他の波数成分同士が混ざらないことを意味し、結局は単一の平面波成分のみを持つという仮定が正しいことなる。
ただし、最初から単色波を仮定する方法は、天下りであり、三階微分以降の非調和項を考えるときにはその導出は全く役に立たなくなる。

2024-07-29

tLeaPで水分子だけのモデルをつくる

AMBER

MDで最も簡単（と思われる）水分子だけのモデルを作りたかったが、ググっても全然出て来なかった。
タンパク質を読み込んで、そこに solvateBOX で水分子を足すことはたくさんやられている。
しかし solvateBOX は "unit" に水分子を加えるため、そもそもの unit を定義する必要がある。
水分子だけを含む unit をパッと作る方法を見つけたので、以下に tLeaP のインプットファイルを記載する。

# tleap.in
source leaprc.water.tip3p
tip3p = sequence { WAT }
#solvateBOX tip3p tip3pBOX 10.0  # All the characters in "TIP3PBOX" should be capital
solvateBOX tip3p TIP3PBOX 10.0
saveAmberParm tip3p tip3p.prmtop tip3p.inpcrd
savePDB tip3p tip3p.pdb
quit

水分子の resname は "WAT" のため、sequence で WAT を１つだけ含む unit を定義すれば、solvateBOX で水分子を追加することができる。

作成した、tleap.in は以下のコマンドで tLeaP に読み込ませられる。

tleap -f tleap.in

2024-02-09

pythonでヒストグラムを作成し、その半値全幅を算出する。

プログラミング数学統計力学

参考サイト
NumPyでヒストグラムを作るnp.histogram関数の使い方 - DeepAge
scipy.signal.peak_widths — SciPy v1.12.0 Manual

import numpy as np
from scipy.signal import find_peak, peak_widths
from matplotlib import pyplot as plt

# data file: data.dat
# data.dat has 2 raw data, and second raw data is to be plotted as a histogram

# read data from data file
data = np.loadtxt('data.dat').T

# check histogram by plotting
plt.hist(data[1])
plt.show()

# get half maximum of full width of the histogram
histogram = np.histgram(data[1])
# histogram[0]: intensity (y-axis)
# histogram[1]: values (x-axis)
peaks = find_peaks(histogram[0])
results_half = peak_widths(histogram[0], peaks, rel_height=0.5)
half_widths = results_half[0] * (histgram[1][1] - histogram[1][0])

# check peak and widths by plotting
plt.plot(histogram[0])
plt.plot(peaks, histogram[0][peaks], "x")
plt.hlines(*results_half[1:])
plt.show()

2024-01-17

sedで、不要な行を消去しながら、範囲を指定して置換し、結果を上書きする

プログラミング

参考サイト
sed で指定した範囲の行を置換する - まくまくsed/awkノート
 sed | テキストの置換処理を得意とするスクリプト言語
 sedでこういう時はどう書く? #Linux - Qiita

text.txt内の1~10行目を無視したい。

sed -e '1,10d' text.txt

無視した後の残された行を、頭から1行目と数え直して、100行目から1000行目の間にある全てのhogeをHOGEに置換したい。

sed -e '1,10d' -e '100,1000 s/hoge/HOGE/g' text.txt

上記の内容を、同じtext.txtに反映させたい。

sed -i -e '1,10d' -e '100,1000 s/hoge/HOGE/g' text.txt

※giにすると、大文字と小文字を区別しなくなる。

#を含む行を消したい。

sed -e '/#/ d' text.txt

行頭に#がある行を消したい。

sed -e '/^#/ d' text.txt

行頭の#を消したい。

sed -e 's/^#//' text.txt

行頭に#がある行のhogeだけ全てHOGEに置換したい。

sed -e '/^#/ s/hoge/HOGE/g' text.txt

- 行頭に#がない行のhogeだけ全てHOGEに置換したい。
>|bash|
sed -e '/^#/! s/hoge/HOGE/g' text.txt

2023-10-23

オイラー・ラグランジュ方程式の導出におけるラグランジュアンの偏微分

力学数学

ラグランジュアンの微分について、佐久間さんの呟きを見つけた。

∂L/∂vに物理的解釈なんか求めなくていいけど、考えずにはいられない物理学徒のために言うと、「関数yを動かすときy,y'は独立ではない。関数を考えるな。架空の位置wと速度vの対(状態)全体のなす空間ではw,vは独立だからwを止めながらvを動かすことが意味をもつ。架空の速度によるLの微分が∂L/∂v」
— 佐久間 (@keisankionwykip) 2021年12月5日

昔の自分の記事で考察したことがあったのを思い出した。
koideforest.hatenadiary.com
改めて見直すと、汎関数微分を全微分で考えようとしていたり、不正確なところもあり、反省した。
その上で、数学と物理の考え方の違いが自分の中で見えた気がする。

ラグランジュアン $L(q(t),\dot q(t),t)$ の考え方に、二つあると見た。

物理：最終的には時間の関数と思う。
数学：独立な２変数（もしくは３変数）を引数に持つ関数と思う。

$q(t)$ と $\dot q(t)$ は時間の関数だから、 $L(q(t),\dot q(t),t)$ に $q(t)$ と $\dot q(t)$ の具体的な形を入れると、結局は時間変数 $t$ だけで表されることになる。
荒っぽく言えば、時間だけの関数を都合の良いように $(q(t),\dot q(t),(t))$ で分けて、それぞれを独立だと思って偏微分すると、それらしいものが得られる。

一方で、本来、 $L$ を関数として考えるならば、定義域をしっかり規定する必要がある。 $L(q,\dot q,(t))$ の形にする時点で、独立な変数の組 $(w,v,\tau)$ の全領域（もしくは物理として必要な領域）で $L(w,v,(\tau))$ が定義されていて、その微分（導関数）も（物理として使う上で）定義されている必要がある。
この導関数は、独立な変数の組 $(w,v,\tau)$ に対する微分操作によって得られる。
物理として $L(q(t),\dot q(t),(t))$ やその導関数を使う場合には、独立な変数の組 $(w,v,\tau)$ の値として $(q(t),\dot q(t),(t))$ を代入したと考えるべきなのだろう。

熱力学の偏微分も、同じようなことが起きていると思われる。

2023-10-20

単純な（熱力学）エントロピーのモデル

単純なエントロピーのモデルを弄ってみる。
参考文献：

清水明、『熱力学の基礎（第２版）』東京大学出版会
H.B. Callen, Thermodynamics and an introduction to thermostatistics, 2nd edition (Wiley, 1985)

エントロピー表示
$\displaystyle S(U,V,N) = (UVN)^{1/3} \, (U \ge 0)$

同次性

$\displaystyle S(\lambda U,\lambda V,\lambda N) = (\lambda^3 UVN)^{1/3} = \lambda (UVN)^{1/3} = \lambda S(U,V, N)$

凸性（上に凸）

$\displaystyle \frac{ \partial^2 S}{\partial U^2} = -\frac{2}{9} \left(\frac{VN}{U^5}\right)^{1/3} \le 0$

エネルギー表示
$\displaystyle U(S,V,N) = \frac{S^3}{VN}$

同次性

$\displaystyle U(\lambda S,\lambda V,\lambda N) = \frac{\lambda^3}{\lambda^2} \frac{S^3}{VN}= \lambda U(S,V, N)$

凸性（下に凸）

$\displaystyle \frac{ \partial^2 U}{\partial S^2} = 6 \frac{S}{VN} \ge 0 \\ \displaystyle \frac{ \partial^2 U}{\partial V^2} = 2 \frac{S}{V^3N} \ge 0$

温度
$\displaystyle \frac{ \partial U}{\partial S} = 3 \frac{S^2}{VN} = \frac{3}{S} U = T \\ \displaystyle U(S, V, N) =\frac{S}{3} T(S, V, N) \\ \displaystyle T(S, V, N) =3\frac{S^2}{VN} \,\, \therefore S = \left(\frac{T}{3} V N\right)^{1/2}$

ヘルムホルツの自由エネルギー
$\displaystyle F(T,V,N) = U(T,V,N) - T S(T,V,N) = -\frac{2}{3} T S(T, V, N) \\ \displaystyle \quad = -\frac{2}{3} \left(\frac{T^3}{3} V N\right)^{1/2}$

凸性（温度に対して上に凸）

$\displaystyle \frac{ \partial^2 F}{\partial T^2} = -\frac{1}{2} \left(\frac{1}{3T} V N\right)^{1/2} \le 0$
温度一定の条件を課せば、凸性は内部エネルギーと一致する。
このことは、凸性がルジャンドル変換で逆になっていても、熱浴下による温度一定条件であれば、平衡状態は $U$ と同様に $F$ を最小することで得られることに対応している。

nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。