速度相関関数とフォノン状態密度

速度自己相関関数（velocity autocorrelation function; VACF）は，次のように書かれる量で，各粒子の異なる時刻との自己相関関数．

$$Z(t)=\frac{⟨v_i(t+t_0)\cdot v_i(t_0)⟩}{⟨v_i(t_0)\cdot v_i(t_0)⟩}=\frac{{\sum }_{t_0}{\sum }_i^N{v}_i(t+t_0)\cdot v_i(t_0)}{{\sum }_{t_0}{\sum }_i^N{v}_i(t_0)\cdot v_i(t_0)}$$

フォノンの状態密度は $Z(t)$ の Fourier 変換として定義される．$Z(t)$ の時間に対する対称性を考慮し，ダンピング関数 $\exp [-(t/{\tau }_{\mathrm{damp}}{)}^2]$ を掛けて，次の形として計算することができる．

$$\tilde{Z}(\omega )=2{\int }_0^{t_{\max }}Z(t)\exp \left[-(t/{\tau }_{\mathrm{damp}}{)}^2\right]\cos (\omega t)\mathrm{d}t,$$

これを離散化して次のように計算する．

$$\tilde{Z}(\omega )=2\sum _{n=0}^{N-1}Z(n\Delta t)\exp \left[-(n\Delta t/{\tau }_{\mathrm{damp}}{)}^2\right]\cos (\omega n\Delta t)\Delta t$$

Phone DOS（もしくは vibrational DOS）を計算する際，次のような振動数領域と時間領域の逆の関係を把握する必要がある．

• 時間刻み → 最大周波数
• 最大時間経過 → 周波数刻み（周波数解像度）

そのため，固体で一般的な周波数上限の 50 THz までの領域のデータを得たいとすると，1/(50 THz) = 20 fs 程度の時間間隔の速度データが必要となる．また，その解像度を 0.1 THz 程度とすると，10 ps 程度の時間区間のデータが必要となる．つまり，この場合だと 500 サンプル程度は必要となる．

上記のシミュレーションをpmdで実行して pmdtraj ファイルに出力されたとすると，napy/nappy/vacf.py を用いて，次のコマンドで VACF および phonon DOS を計算することができる．

python ~/src/nap/nappy/vacf.py --normalize --sigma 5 -t 20.0 pmdtraj

• -t 20： データ間の時間刻み（fs）．
• --sigma： 小さな系で得られた DOS はガタガタしていることが多いため，ここで指定した周波数点の幅のガウシアンでぼかす．
• --normalize： phonon DOS を全体のDOSの積分が１となるように規格化する．